/**************************************************** * Reduction of ternary forms * * Michael Stoll, started 2002-02-21 ****************************************************/ /**************************************************** We associate to a nonsingular ternary form C(x,y,z) a positive definite quadratic form Q(x,y,z). This is done as follows. When C is a nonsingular cubic: Let H(x,y,z) = det(matrix of second partial derivatives of C) be the Hessian of C. Then the intersection of C and H consists of 9 distinct points, the flexes of C. Three of them are real, the others come in three complex conjugate pairs. There are twelve lines each containing three of teh flexes, coming in four triples of lines that do not meet in a flex. One of these triples has all three lines real, call them L_11, L_12, L_13. Another one has one line real, call it L_21, and two complex conjugate lines, call them L_22 and L_23. Then our form Q spans the one-dimensional intersection of the spaces spanned by L_11^2, L_12^2 and L_13^2, and by L_21^2 and L_22 L_23, respectively. This Q is positive definite and (as a Hermitian form) a covariant of C (up to scaling); in fact the unique (real) Hermitian covariant of C. If C is a higher degree form, we still compute the flexes (3 d (d-2) of them, where d is the degree of C). Now we use them to find a covariant positive definite Hermitian (or real quadratic) form Q by some optimisation procedure. Q defines a lattice which we can Minkowsi- or LLL-reduce. Applying the corresponding transformations to C, we obtain a reduced curve. I.e., C is (LLL-)reduced if Q is reduced if the lattice defined by Q is Minkowski-reduced (LLL-reduced). ****************************************************/ declare verbose TernaryReduce, 3; declare verbose Minimise, 3; intrinsic PolMat(pol::RngMPolElt, mat::Mtrx) -> RngMPolElt {pol^mat; workaround for bug in Magma.} P := Parent(pol); n := Rank(P); require n eq NumberOfRows(mat) and n eq NumberOfColumns(mat): "mat must be a square matrix of order equal to the number of variables in pol."; return Evaluate(pol, [&+[mat[i,j]*P.j : j in [1..n]] : i in [1..n]]); end intrinsic; // Generate random matrices in SL(n, Z) intrinsic RandomSLnZ(n::RngIntElt, k::RngIntElt, l::RngIntElt) -> AlgMatElt {A random element in SL(n,Z). k and l have an influence on the size of the entries.} eltmat := func; return &*[eltmat(Random(1,n), Random(1,n), Random(-k,k)) : j in [1..l]]; end intrinsic; // Generate random forms intrinsic RandomForm(P::RngMPol, deg::RngIntElt, k::RngIntElt) -> RngMPolElt {Generates a random form of degree deg in P. k bounds the coefficients.} return &+[Random(-k,k)*m : m in MonomialsOfDegree(P, deg)]; end intrinsic; intrinsic ComplexVariety(C::RngMPolElt, H::RngMPolElt) -> SeqEnum {Computes the complex points of the zero-dimensional projective variety given by the two ternary forms C and H over Q.} P := Parent(C); Q := CoefficientRing(P); require Type(Q) eq FldRat: "Base field must be the rationals."; require P eq Parent(H): "Both forms must have the same parent."; require Rank(P) eq 3: "Forms must be in three variables."; dC := TotalDegree(C); dH := TotalDegree(H); require forall{m : m in Monomials(C) | TotalDegree(m) eq dC} and forall{m : m in Monomials(H) | TotalDegree(m) eq dH}: "Polynomials must be forms (i.e., homogeneous)."; require TotalDegree(GCD(C, H)) eq 0: "Curves have a component in common."; PC := PolynomialRing(ComplexField()); Proj := ProjectiveSpace(Q, 2); PP := CoordinateRing(Proj); CC := Curve(Proj, PP!C); CH := Curve(Proj, PP!H); if Evaluate(C, [1,0,0]) eq 0 and Evaluate(H, [1,0,0]) eq 0 then pt := Proj![1,0,0]; points := [<[1,0,0], IntersectionNumber(CC, CH, pt)>]; else points := []; end if; PQ := PolynomialRing(Q); res := GCD(Evaluate(C, [PQ.1,1,0]), Evaluate(H, [PQ.1,1,0])); fact := Factorisation(res); fact := [ where PrK := BaseChange(Proj, K) where r := Evaluate(res, K.1) eq 0 select K.1 else Roots(res, K)[1,1] where K := ext : f in fact]; points cat:= &cat[[<[r[1], 1, 0], f[2]> : r in Roots(f[1], ComplexField())] : f in fact]; res := Evaluate(Resultant(Evaluate(C, [P.1,P.2,1]), Evaluate(H, [P.1,P.2,1]), P.2), [PQ.1,0,0]); fact := Factorisation(res); for f in fact do f1 := f[1]; if Degree(f1) eq 1 then r := -Coefficient(f1,0); res1 := GCD(Evaluate(C, [r, PQ.1, 1]), Evaluate(H, [r, PQ.1, 1])); fact1 := SquarefreeFactorization(res1); for ff in fact1 do points cat:= [<[r, rs[1], 1], ff[2]> : rs in Roots(ff[1], ComplexField())]; end for; elif Degree(f1) le 6 or f[2] ge 2 then // Do computations exactly as far as possible // to get muliplicities right. K := ext; if Evaluate(f1, K.1) eq 0 then r := K.1; else r := Roots(f1, K)[1,1]; end if; PK := PolynomialRing(K); res1 := GCD(Evaluate(C, [r, PK.1, 1]), Evaluate(H, [r, PK.1, 1])); fact1 := SquarefreeFactorization(res1); for rp in Roots(f1, ComplexField()) do h := hom ComplexField() | rp[1]>; for ff in fact1 do points cat:= [<[rp[1], r[1], 1], ff[2]> : r in Roots(PC![h(c) : c in Coefficients(ff[1])])]; end for; end for; else // Try to get by without working in high degree number fields. // Result may be less accurate, though. _, prec := HasAttribute(FldPr, "Precision"); oldprec := prec; repeat ptsloc := []; flag := true; for r in [r[1] : r in Roots(f1, ComplexField())] do pol1 := Evaluate(C, [r, PC.1, 1]); pol2 := Evaluate(H, [r, PC.1, 1]); // We know that pol1 and pol2 have a unique common root. // Do the euclidean algorithm until we are left with a linear factor. pol1 *:= 1/LeadingCoefficient(pol1); pol2 *:= 1/LeadingCoefficient(pol2); p1 := pol1; p2 := pol2; while Degree(pol2) gt 1 do temp := pol2 mod pol1; if temp eq 0 then break; end if; pol2 := pol1; if Degree(pol2) gt 1 then pol1 := 1/LeadingCoefficient(temp)*temp; end if; end while; r2 := -Coefficient(pol2, Degree(pol2)-1)/Degree(pol2); if Norm(Evaluate(p1, r2)) ge 1.0e-6 or Norm(Evaluate(p2, r2)) ge 1.0e-6 then prec +:= 50; AssertAttribute(FldPr, "Precision", prec); vprintf TernaryReduce, 2: "ComplexVariety: Increasing precision to %o\n", prec; flag := false; break; end if; Append(~ptsloc, <[ComplexField()|r, r2, 1], 1>); end for; until flag; AssertAttribute(FldPr, "Precision", oldprec); points cat:= ptsloc; end if; end for; return points; end intrinsic; // Find real and complex flexes of a ternary form over Q intrinsic RealAndComplexFlexes(C::RngMPolElt) -> SeqEnum, SeqEnum {Finds the real and complex flexes of a ternary form over Q. First value is the sequence of real flexes, second value is a sequence containing one complex flex out of each conjugate pair.} P := Parent(C); deg := TotalDegree(C); require Rank(P) eq 3 and forall{m : m in Monomials(C) | TotalDegree(m) eq deg}: "C must be a ternary form."; vprintf TernaryReduce, 1: "Flexes: C = %o\n", C; // compute Hessian H := Determinant(Matrix([[Derivative(Derivative(C,i),j) : j in [1..3]] : i in [1..3]])); vprintf TernaryReduce, 2: " Hessian H = %o\n", H; flexes := &cat[[v[1] : i in [1..v[2]]] : v in ComplexVariety(C, H)]; require #flexes eq 3*deg*(deg-2): "Did not find the correct number of flexes."; flr := [fl : fl in flexes | IsReal(fl[1]) and IsReal(fl[2])]; flc0 := [fl : fl in flexes | not IsReal(fl[1]) or not IsReal(fl[2])]; // find representatives of conjugate pairs flc := []; for i := 1 to ExactQuotient(#flc0, 2) do new := flc0[#flc0]; cnew := [ComplexConjugate(new[i]) : i in [1..3]]; Append(~flc, new); min, pos := Minimum([&+[Norm(cnew[i]-flc0[j][i]) : i in [1..3]] : j in [1..#flc0-1]]); assert min le 1.0e-5*&+[Norm(new[i]) : i in [1..3]]; // assert pos ne #flc0; Remove(~flc0, pos); Prune(~flc0); end for; vprintf TernaryReduce, 3: "Flexes: %o\n", flexes; return flr, flc; end intrinsic; // Write my own Kernel function, since I got a ``division by zero in gdiv'' function MyKernel(mat) // mat = [[a,b], [c,d], [e,f]] // find linear combination of rows that vanishes // so operate on columns V := KSpace(CoefficientRing(mat), 3); // find largest entry in matrix max, pos := Max([Abs(c) : c in Eltseq(mat)]); if IsEven(pos) then // 2nd col SwapColumns(~mat, 1, 2); pos -:= 1; end if; r := ExactQuotient(pos+1, 2); // row number mat *:= 1/mat[r,1]; AddColumn(~mat, -mat[r,2], 1, 2); // rth col is (1 0) // find largest entry in 2nd col max, pos := Max([Abs(mat[i,2]) : i in [1..3]]); if pos eq r or max eq 0 then return Eltseq([V.i - mat[i,1]*V.r : i in [1..3] | i ne r][1]); end if; MultiplyColumn(~mat, 1/mat[pos,2], 2); AddColumn(~mat, -mat[pos,1], 2, 1); // pos'th col is (0 1) t := 6-r-pos; // the third element of {1,2,3} = {r,pos,t} return Eltseq(V.t - mat[t,1]*V.r - mat[t,2]*V.pos); end function; // Given the flexes of a ternary cubic, find the relevant lines function LinesFromFlexes(freal, fcomplex) // Given sequences containing the three real flexes and three non-conjugate // complex flexes of a ternary cubic, respectively, finds the // line through the 3 real flexes, its complementary factor, and // the three real lines through a real and two complex flexes. P3 := PolynomialRing(RealField(), 3); normalise := function(l) _, pos := Maximum([Norm(a) : a in l]); norm := Sqrt(&+[Norm(a) : a in l]) * l[pos]/Modulus(l[pos]); return [a/norm : a in l]; end function; line := function(pt1, pt2) l := MyKernel(Matrix([[pt1[i], pt2[i]] : i in [1..3]])); return normalise(l); end function; lines := [line(freal[1],freal[2]), line(freal[2],freal[3]), line(freal[3],freal[1])]; liner := normalise([&+[l[i] : l in lines] : i in [1..3]]); vprintf TernaryReduce, 3: " Real line: %o\n", liner; linec := [[Real(a) : a in line(fl, [ComplexConjugate(c) : c in fl])] : fl in fcomplex]; vprintf TernaryReduce, 3: " Real triple of lines:\n%o\n", linec; linetest1 := line(fcomplex[1], fcomplex[2]); linetest2 := line(fcomplex[1], [ComplexConjugate(c) : c in fcomplex[2]]); vprintf TernaryReduce, 3: " Possible complex lines:\n %o\n %o\n", linetest1, linetest2; m1 := Min(Norm(&+[linetest1[i]*fcomplex[3,i] : i in [1..3]]), Norm(&+[linetest1[i]*ComplexConjugate(fcomplex[3,i]) : i in [1..3]])); m2 := Min(Norm(&+[linetest2[i]*fcomplex[3,i] : i in [1..3]]), Norm(&+[linetest2[i]*ComplexConjugate(fcomplex[3,i]) : i in [1..3]])); error if Min(m1, m2) ge 1.0e-10, "Lines do not fit!"; vprintf TernaryReduce, 3: " Values of two possible lines at third point: %o, %o\n", m1, m2; linet := m1 lt m2 select linetest1 else linetest2; compl := Norm(linet[1])*P3.1^2 + Norm(linet[2])*P3.2^2 + Norm(linet[3])*P3.3^2 + 2*Real(linet[1]*ComplexConjugate(linet[2]))*P3.1*P3.2 + 2*Real(linet[1]*ComplexConjugate(linet[3]))*P3.1*P3.3 + 2*Real(linet[2]*ComplexConjugate(linet[3]))*P3.2*P3.3; vprintf TernaryReduce, 3: " Complementary factor:\n%o\n", compl; return [&+[liner[i]*P3.i : i in [1..3]], compl] cat [&+[l[i]*P3.i : i in [1..3]] : l in linec]; end function; function orthonorm(basis) for i := 1 to #basis do s := [&+[basis[j,k]*basis[i,k] : k in [1..#basis[i]]] : j in [1..i-1]]; if i gt 1 then basis[i] := [basis[i,k] - &+[s[j]*basis[j,k] : j in [1..#s]] : k in [1..#basis[i]]]; end if; n := Sqrt(&+[a^2 : a in basis[i]]); basis[i] := [a/n : a in basis[i]]; end for; return basis; end function; // Now find the form from the lines function CovariantFromLines(lines) // Finds the covariant positive definite Hermitian form from the // given lines through the flexes of a ternary cubic. P := Universe(lines); // look for element in space that is closest to // mons := MonomialsOfDegree(P, 2); w1 := [ MonomialCoefficient(lines[1]^2, mon) : mon in mons ]; w2 := [ MonomialCoefficient(lines[2], mon) : mon in mons ]; // orthonormalise vprintf TernaryReduce, 3: " First basis:\n %o\n %o\n", w1, w2; w := orthonorm([w1, w2]); vprintf TernaryReduce, 3: " orthonormalised:\n %o\n %o\n", w[1], w[2]; // basis of other space v1 := [ MonomialCoefficient(lines[3]^2, mon) : mon in mons ]; v2 := [ MonomialCoefficient(lines[4]^2, mon) : mon in mons ]; v3 := [ MonomialCoefficient(lines[5]^2, mon) : mon in mons ]; // orthonormalise vprintf TernaryReduce, 3: " Second basis:\n %o\n %o\n %o\n", v1, v2, v3; v := orthonorm([v1, v2, v3]); vprintf TernaryReduce, 3: " orthonormalised:\n %o\n %o\n %o\n", v[1], v[2], v[3]; // look at cos(phi)*w[1] + sin(phi)*w[2] such that distance to // is minimal. Leads to scmat := [[&+[wi[k]*vj[k] : k in [1..#mons]] : wi in w] : vj in v]; vprintf TernaryReduce, 3: " Sclar products matrix:\n%o\n", scmat; phi := 1/2*Arctan2(&+[a[1]*a[2] : a in scmat], 1/2*&+[a[1]^2 - a[2]^2 : a in scmat]); cphi := Cos(phi); sphi := Sin(phi); qs1 := [cphi*w[1,k] + sphi*w[2,k] : k in [1..#mons]]; qs2 := [-sphi*w[1,k] + cphi*w[2,k] : k in [1..#mons]]; ls1 := [cphi*a[1] + sphi*a[2] : a in scmat]; ls2 := [-sphi*a[1] + cphi*a[2] : a in scmat]; dist1 := 1 + &+[l^2 : l in ls1] - 2*&+[ls1[i]*(cphi*scmat[i,1] + sphi*scmat[i,2]) : i in [1..#ls1]]; dist2 := 1 + &+[l^2 : l in ls2] - 2*&+[ls2[i]*(-sphi*scmat[i,1] + cphi*scmat[i,2]) : i in [1..#ls2]]; if dist1 lt dist2 then qs := qs1; ls := ls1; dist := dist1; else qs := qs2; ls := ls2; dist := dist2; end if; vprintf TernaryReduce, 3: " Nearest vecor in first space:\n %o\n", qs; vv := [&+[ls[i]*v[i,k] : i in [1..#ls]] : k in [1..#mons]]; vprintf TernaryReduce, 3: " Nearest vector in second space:\n %o\n", vv; Q := &+[qs[i]*mons[i] : i in [1..#mons]]; if MonomialCoefficient(Q, Parent(Q).1^2) lt 0 then Q := -Q; end if; Qmat := Matrix([[(i eq j select 1 else 1/2)*MonomialCoefficient(Q, P.i*P.j) : j in [1..3]] : i in [1..3]]); error if not IsPositiveDefinite(Qmat), "Covariant matrix is not positive definite -- try with higher precision."; Qmat *:= Determinant(Qmat)^(-1/3); return Qmat, dist; end function; intrinsic Covariant(C::RngMPolElt : Al := "Default", Eps := 1.0e-20) -> RngMPolElt {Computes the positive definite Hermitian covariant of a ternary cubic as a real quadratic form.} require Rank(Parent(C)) eq 3: "C must be a ternary form."; deg := TotalDegree(C); require forall{m : m in Monomials(C) | TotalDegree(m) eq deg}: "Polynomial must be homogeneous."; require deg ge 3: "Degree must be at least 3."; require Al in {"Default", "SteepestDescent", "Newton"}: "Al must be one of \"Default\", \"SteepestDescent\" and \"Newton\""; if deg eq 3 and Al eq "Default" then flr, flc := RealAndComplexFlexes(C); return CovariantFromLines(LinesFromFlexes(flr, flc)); else flr, flc := RealAndComplexFlexes(C); return Covariant(flr, flc : Al := Al, Eps := Eps); end if; end intrinsic; /////////////////////////////////////////////////////////////////////// // Ad-hoc reduction functions /////////////////////////////////////////////////////////////////////// function siz(pol) return Max([Abs(c) : c in Coefficients(pol)]); end function; function redstep(cubic) P := Parent(cubic); z := [P.1, P.2, P.3]; trans := IdentityMatrix(Integers(), 3); for i := 1 to 3 do mc := MonomialCoefficient(cubic, z[i]^3); if mc ne 0 then as := [j eq i select 1 else -Round(MonomialCoefficient(cubic, z[i]^2*z[j])/(3*mc)) : j in [1..3]]; tr := Matrix([ii eq i select as else [j eq ii select 1 else 0 : j in [1..3]] : ii in [1..3]]); cubic := PolMat(cubic, tr); trans := trans*tr; end if; end for; return cubic, trans; end function; function redrep(C) trans := IdentityMatrix(Integers(), 3); s0 := siz(C); s1 := siz(C); C1 := C; tr1 := trans; repeat s := s1; C := C1; trans := tr1; C1, tr := redstep(C); tr1 := trans*tr; s1 := siz(C1); until s1 ge s; return C, trans; end function; function em(i,j,k) // elementary matrix return IdentityMatrix(Integers(),3) + Matrix(3, [(ii eq i and jj eq j) select k else 0 : ii,jj in [1..3]]); end function; // set up a list of (monoid) generators gs1 := [em(i,j,1) : i,j in [1..3] | i ne j]; gs2 := [em(i,j,-1) : i,j in [1..3] | i ne j]; gens1 := gs1 cat gs2; gens2 := gens1 cat [gs1[i]*gs1[j] : j in [i..#gs1], i in [1..#gs1]] cat [gs1[i]*gs2[j] : i,j in [1..#gs1] | i ne j] cat [gs2[i]*gs2[j] : j in [i..#gs2], i in [1..#gs2]]; function red0(pol, gens) pols := [PolMat(pol, g) : g in gens]; _, pos := Min([siz(mm) : mm in pols]); return pols[pos], gens[pos]; end function; function red0rep(C) trans := IdentityMatrix(Integers(), 3); s1 := siz(C); C1 := C; tr1 := trans; dig := Ceiling(Log(10, s1)); repeat s := s1; C := C1; trans := tr1; C1, tr := red0(C, gens2); tr1 := trans*tr; s1 := siz(C1); dig1 := Ceiling(Log(10, s1)); if dig1 lt dig then vprintf TernaryReduce, 2: " red0rep: size now = %o digits\n", dig1; dig := dig1; end if; until s1 ge s; return C, trans; end function; function AdHocRedCubic(C) trans := IdentityMatrix(Integers(), 3); s1 := siz(C); C1 := C; tr1 := trans; dig := Ceiling(Log(10, s1)); count := 100; repeat count -:= 1; s := s1; C := C1; trans := tr1; C1, tr := redrep(C); tr1 := trans*tr; C1, tr := red0(C1, gens2); tr1 := tr1*tr; s1 := siz(C1); dig1 := Ceiling(Log(10, s1)); vprintf TernaryReduce, 3: " Ad-hoc: gain = %o\n", s-s1; if dig1 lt dig then vprintf TernaryReduce, 2: " Ad-hoc cubic red: size now = %o digits\n", dig1; dig := dig1; end if; until s1 ge s or count le 0; return C, trans; end function; //////////////////////////////////////////////////////////////////////// intrinsic Reduce(C::RngMPolElt : Al := "Default", Steps := 3, AdHoc := true) -> RngMPolElt, AlgMatElt {Computes an LLL-reduced form in the SL_3(Z)-equivalence class of C.} P := Parent(C); require Rank(P) eq 3: "C must be a ternary form."; deg := TotalDegree(C); require forall{m : m in Monomials(C) | TotalDegree(m) eq deg}: "Polynomial must be homogeneous."; require deg ge 3: "Degree must be at least 3."; C0 := C; vprintf TernaryReduce, 1: "Reduce: C = %o\n", C; P1 := PolynomialRing(BaseRing(P)); PR := PolynomialRing(RealField()); trans := IdentityMatrix(Integers(), 3); // do some preliminary steps if AdHoc then repeat if deg eq 3 then vprintf TernaryReduce, 2: " Preliminary ad-hoc reduction for cubics...\n"; C1, tr := AdHocRedCubic(C); trans := trans*tr; else vprintf TernaryReduce, 2: " Preliminary ad-hoc reduction for general curves...\n"; C1, tr := red0rep(C); trans := trans*tr; end if; ff := C1 eq C; C := C1; vprintf TernaryReduce, 2: " Ad-hoc reduction...\n"; deg1 := 3*deg*(deg-2); count := 10; repeat s := siz(C); H := TernaryHessian(C); vprintf TernaryReduce, 2: " Computing resultants...\n"; s1 := siz(C); BI1 := Resultant(C, H, P.3); BI2 := Resultant(C, H, P.2); BI3 := Resultant(C, H, P.1); pol12 := Evaluate(BI1, [P1.1,1,0]); pol21 := Evaluate(BI1, [1,P1.1,0]); pol13 := Evaluate(BI2, [P1.1,0,1]); pol31 := Evaluate(BI2, [1,0,P1.1]); pol23 := Evaluate(BI3, [0,P1.1,1]); pol32 := Evaluate(BI3, [0,1,P1.1]); pols := [pol12,pol21,pol13,pol31,pol23,pol32]; pols := [pol/LeadingCoefficient(pol) : pol in pols]; // see if complex roots are close together eps := 1/1000; rtsav := [Degree(pol) lt deg1 select 0 else -Coefficient(pol,deg1-1)/deg1 : pol in pols]; rtsr := [Round(r) : r in rtsav]; flags := [rtsr[i] ne 0 and forall{j : j in [0..deg1-1] | Abs(Coefficient(pp,j)) le Binomial(deg1,j)*eps^(deg1-j)} where pp := Evaluate(pols[i], P1.1+rtsav[i]) : i in [1..6]]; mats := []; if flags[1] and flags[6] then // affine y = 1 Append(~mats, Matrix([[1,rtsr[1],0], [0,1,0], [0,rtsr[6],1]])); end if; if flags[2] and flags[4] then // affine x = 1 Append(~mats, Matrix([[1,0,0], [rtsr[2],1,0], [rtsr[4],0,1]])); end if; if flags[3] and flags[5] then // affine z =1 Append(~mats, Matrix([[1,0,rtsr[3]], [0,1,rtsr[5]], [0,0,1]])); end if; if not IsEmpty(mats) then _, pos := Max([&+[Abs(c) : c in Eltseq(m)] : m in mats]); tr := mats[pos]; elif &or flags then if flags[1] then Append(~mats, Matrix([[1,rtsr[1],0], [0,1,0], [0,0,1]])); end if; if flags[2] then Append(~mats, Matrix([[1,0,0], [rtsr[2],1,0], [0,0,1]])); end if; if flags[3] then Append(~mats, Matrix([[1,0,rtsr[3]], [0,1,0], [0,0,1]])); end if; if flags[4] then Append(~mats, Matrix([[1,0,0], [0,1,0], [rtsr[4],0,1]])); end if; if flags[5] then Append(~mats, Matrix([[1,0,0], [0,1,rtsr[5]], [0,0,1]])); end if; if flags[6] then Append(~mats, Matrix([[1,0,0], [0,1,0], [0,rtsr[6],1]])); end if; _, pos := Max([&+[Abs(c) : c in Eltseq(m)] : m in mats]); tr := mats[pos]; else break; end if; ff := false; C := PolMat(C, tr); trans := trans*tr; s1 := siz(C); vprintf TernaryReduce, 2: " Ad-hoc red: Size now %o digits\n", Ceiling(Log(10, s1)); if s1 ge s then count -:= 1; else s := s1; count := 5; end if; until count eq 0; until ff; end if; //AdHoc assert C eq PolMat(C0, trans); vprintf TernaryReduce, 2: " After ad-hoc reduction, curve is\n%o\n", C; // now do the serious work for eps in [1.0e-5, 1.0e-10, 1.0e-20][1..Max(1,Min(Steps,3))] do vprintf TernaryReduce, 2: " Reduce: Using eps = %o\n", eps; cov := Covariant(C : Al := Al, Eps := eps); _, tr := LLLGram(cov); tr := Transpose(tr); vprintf TernaryReduce, 2: " Transformation:\n%o\n", tr; C := PolMat(C, tr); vprintf TernaryReduce, 2: " New C = %o\n", C; trans := trans*tr; end for; // reduce completely mins1 := []; c := 16; repeat mins := mins1; repeat c *:= 2; M := Matrix(3,[Round(c*e) : e in Eltseq(cov)]); until IsPositiveDefinite(M); _, mins1 := SuccessiveMinima(LatticeWithGram(M)); mins1 := [Eltseq(v) : v in mins1]; until mins1 eq mins; tr := Transpose(Matrix(mins)); vprintf TernaryReduce, 2: " Final transformation:\n%o\n", tr; C := PolMat(C, tr); vprintf TernaryReduce, 2: " Final C = %o\n", C; return C, trans*tr; end intrinsic; // Now try general ternary forms. // Use the flexes to find a covariant point geometrically // Substitute for missing Exp(mat) function SimpleMatExp(mat) n := NumberOfRows(mat); e := Parent(mat)!IdentityMatrix(RealField(), n); mat -:= 1/n*Trace(mat)*e; norm := Sqrt(&+[Norm(c) : c in Eltseq(mat)]); b := Ceiling(10*norm); mat *:= 1/b; m := mat; for j := 1 to 10 do e +:= m; m *:= 1/(j+1)*mat; end for; e +:= m; e := e^b; e *:= (Determinant(e))^(-1/n); return e; end function; intrinsic Covariant(ptseqr::SeqEnum, ptseqc::SeqEnum : Steps := 0, Eps := 1.0e-20, Al := "Newton") -> AlgMatElt, FldPrElt {Given two sequences of real, resp. complex, points, represented as n-tuples of coordinates, not all equal to zero, finds the covariant positice definite Hermitian form as a symmetric real matrix. The collection of points (real + complex + conjugate complex) must be balanced. (I.e., less than d/n of the points may be in a d-dimensional subspace for d = 1, ..., n-1. This is not checked; if not satisfied, teh algorithm will not terminate.)} s := {#pt : pt in ptseqr cat ptseqc}; require #s eq 1: "All points must have the same number of coordinates."; n := Rep(s); nr := #ptseqr; nc := #ptseqc; npts := nr + 2*nc; require nr + 2*nc gt n: "Need at least %o points.", n+1; require Al in {"Newton", "SteepestDescent", "Default"}: "Al must be one of \"Default\", \"Newton\" and \"SteepestDescent\""; vprintf TernaryReduce, 1: "Covariant of points:\n Real: %o\n Complex: %o\n", ptseqr, ptseqc; R := RealField(); // work at I trans := IdentityMatrix(R, n); M := Parent(trans); theta := func; norm := func; grad := func; matexp := func; newseqs := func; // Do a first step if points are very close together ptseq := ptseqr cat ptseqc; dth := 0; if forall{pt : pt in ptseq | pt[3] eq 1} then mid := [(&+[pt[i] : pt in ptseqr] + 2*&+[Real(pt[i]) : pt in ptseqc])/npts : i in [1,2]]; rad := Max([Norm(pt[1]-mid[1]) + Norm(pt[2]-mid[2]) : pt in ptseq]); if rad lt 0.1 then // shift mid to (0,0) and scale to make rad = 1 vprintf TernaryReduce, 2: " Points are very close to (%o, %o)\n", mid[1], mid[2]; mrad := (&+[Norm(pt[1]-mid[1]) + Norm(pt[2]-mid[2]) : pt in ptseqr] + 2*&+[Norm(pt[1]-mid[1]) + Norm(pt[2]-mid[2]) : pt in ptseqc]) / npts; r := mrad^(-1/6); trans := Matrix(R, [[1,0,0], [0,1,0], [-mid[1],-mid[2],1]]) * DiagonalMatrix([r,r,r^-2]); vprintf TernaryReduce, 3: " Transforming by\n%o\n", trans; ptseqr, ptseqc := newseqs(ptseqr, ptseqc, trans); ptseqr := [[r^2*e : e in pt] : pt in ptseqr]; ptseqc := [[r^2*e : e in pt] : pt in ptseqc]; dth := npts*4*Log(r); vprintf TernaryReduce, 3: " Points are now\n%o\n%o\n", ptseqr, ptseqc; end if; end if; prec := Min([pr : c in &cat[[pt[1],pt[2],pt[3]] : pt in ptseqr cat ptseqc] | pr ge 0 where pr := Precision(c)]); _, oldprec := HasAttribute(FldPr, "Precision"); Eps := Max(Eps, Exp(-1.5*Log(10)*prec)); prec := Max(oldprec, prec); AssertAttribute(FldPr, "Precision", prec); vprintf TernaryReduce, 2: " Computing with precision = %o\n", prec; // Set up machinery for Newton steps eltmat := func< i,j | Matrix(R, n, [ii eq i and jj eq j select 1/Sqrt(2) else 0 : jj, ii in [1..n]]) >; tr := Orthonormalize(Matrix(R, n-1, [i eq j select 2 else 1 : i,j in [2..n]])); b1mat := tr^-1*Matrix(R, [[i eq 1 select 1 else i eq j select -1 else 0 : i in [1..n]] : j in [2..n]]); basis := [DiagonalMatrix([b1mat[i,j] : j in [1..n]]) : i in [1..n-1]] cat [eltmat(i,j) + eltmat(j,i) : j in [i+1..n], i in [1..n]]; // basis is an orthonormal basis of the tangent space at I ipr := func; ipc := func; function newton(ptseqr, ptseqc) // 1st derivative der1 := [&+[ipr(pt, b)/norm(pt) : pt in ptseqr] + 2*&+[ipc(pt, b)/norm(pt) : pt in ptseqc] : b in basis]; // 2nd derivative der2 := Matrix([[&+[ipr(pt,b1*b2)/norm(pt) - (ipr(pt,b1)*ipr(pt,b2))/norm(pt)^2 : pt in ptseqr] + 2*&+[ipc(pt,b1*b2)/norm(pt) - (ipc(pt,b1)*ipc(pt,b2))/norm(pt)^2 : pt in ptseqc] : b2 in basis] : b1 in basis]); // Newton step step := -1/2*Vector(der1)*der2^-1; gr := &+[der1[i]*basis[i] : i in [1..#basis]]; // gradient return gr, &+[step[i]*basis[i] : i in [1..#basis]]; end function; // Try Newton steps first, do steepest descent if that doesn't reduce theta. c := Exp(0); // get 1.0 to current precision oldtheta := theta(ptseqr, ptseqc); vprintf TernaryReduce, 2: " Initial theta = %o\n", oldtheta-dth; repeat // Newton first flag := false; if Al eq "Newton" then gr, step := newton(ptseqr, ptseqc); gr := -1/2*(gr + Transpose(gr)); vprintf TernaryReduce, 2: " Gradient = \n%o\n", gr; step := 1/2*(step + Transpose(step)); vprintf TernaryReduce, 3: " Newton step =\n%o\n", step; if &+[Norm(c) : c in Eltseq(step)] gt 10 then vprintf TernaryReduce, 3: " Step too large -- try steepest descent.\n"; flag := true; else mat := matexp(step); mat := 1/2*(mat + Transpose(mat)); vprintf TernaryReduce, 3: " Exponentiated =\n%o\n", mat; psr1, psc1 := newseqs(ptseqr, ptseqc, mat); th1 := theta(psr1, psc1); if th1 lt oldtheta then vprintf TernaryReduce, 3: " Newton step successful.\n"; else vprintf TernaryReduce, 3: " Newton step unsuccessful -- try steepest descent.\n"; flag := true; end if; end if; else gr := -grad(ptseqr, ptseqc); vprintf TernaryReduce, 2: " Gradient = \n%o\n", gr; flag := true; end if; if flag then mat := matexp(c*gr); mat := 1/2*(mat + Transpose(mat)); psr1, psc1 := newseqs(ptseqr, ptseqc, mat); th1 := theta(psr1, psc1); while th1 gt oldtheta do c /:= 2; vprintf TernaryReduce, 3: " th1 > th0 : reduce c to %o\n", c; mat := matexp(c*gr); mat := 1/2*(mat + Transpose(mat)); psr1, psc1 := newseqs(ptseqr, ptseqc, mat); th1 := theta(psr1, psc1); end while; mat2 := matexp(c/2*gr); mat2 := 1/2*(mat2 + Transpose(mat2)); psr2, psc2 := newseqs(ptseqr, ptseqc, mat2); th2 := theta(psr2, psc2); if th2 lt th1 then th1 := th2; c /:= 2; psr1 := psr2; psc1 := psc2; mat := mat2; else mat2 := matexp(c*2*gr); mat2 := 1/2*(mat2 + Transpose(mat2)); psr2, psc2 := newseqs(ptseqr, ptseqc, mat2); th2 := theta(psr2, psc2); if th2 lt th1 then th1 := th2; c *:= 2; psr1 := psr2; psc1 := psc2; mat := mat2; end if; end if; vprintf TernaryReduce, 3: " c now = %o\n", c; end if; oldtheta := th1; vprintf TernaryReduce, 2: " Theta now = %o\n", th1 - dth; ptseqr := psr1; ptseqc := psc1; vprintf TernaryReduce, 3: " Transform points by:\n%o\n", mat; trans := trans*mat; vprintf TernaryReduce, 3: " New total transformation:\n%o\n", trans; until &+[c^2 : c in Eltseq(gr)] lt Eps or c^2 lt Eps; AssertAttribute(FldPr, "Precision", oldprec); herm := trans*Transpose(trans); vprintf TernaryReduce, 2: "Covariant is\n%o\n", herm; vprintf TernaryReduce, 2: "Theta = %o\n", oldtheta - dth; return herm, oldtheta - dth; end intrinsic; ////////////////////////////////////////////////////////////////////// // Ueberschiebung(F,G,k) // = c*det([[d/dx1 d/dx2] [d/dy1 d/dy2]])^k // F(x1,y1) G(x2,y2)|(x1=x2=x, y1=y2=y) // Can generalise this to three ternary forms: // TernaryUeberschiebung(F_1,F_2,F_3, k) // = c*det((d/dx_ij))^k F_1(x_11, x_12, x_13) F_2() F_3()|(x_ij = x_j) intrinsic TernaryUeberschiebung(F1::RngMPolElt, F2::RngMPolElt, F3::RngMPolElt, k::RngIntElt) -> RngMPolElt {Computes the k-th Ueberschiebung of three ternary forms.} P3 := Parent(F1); require Parent(F2) eq P3 and Parent(F3) eq P3: "TernaryUeberschiebung(",F1,", ",F2,", ",F3,", ",k,"): ", "Forms must live in the same polynomial ring."; require Rank(P3) eq 3: "TernaryUeberschiebung(",F1,", ",F2,", ",F3,", ",k,"): ", "Forms must be ternary."; x0 := [P3.1, P3.2, P3.3]; n1 := TotalDegree(F1); n2 := TotalDegree(F2); n3 := TotalDegree(F3); require forall{m : m in Monomials(F1) | TotalDegree(m) eq n1} and forall{m : m in Monomials(F2) | TotalDegree(m) eq n2} and forall{m : m in Monomials(F3) | TotalDegree(m) eq n3}: "TernaryUeberschiebung(",F1,", ",F2,", ",F3,", ",k,"): ", "Polynomials must be homogeneous."; if n1 lt k or n2 lt k or n3 lt k then return P3!0; end if; P9<[x]> := PolynomialRing(CoefficientRing(P3), 9); td := func; op := func; res := Evaluate(F1, x[1..3])*Evaluate(F2, x[4..6])*Evaluate(F3, x[7..9]); for i := 1 to k do res := op(res); end for; coef := 1/(&*[n : n in [n1-k+1..n1]] * &*[n : n in [n2-k+1..n2]] * &*[n : n in [n3-k+1..n3]]); return (CoefficientRing(P3)!coef)*Evaluate(res, x0 cat x0 cat x0); end intrinsic; intrinsic TernaryHessian(F::RngMPolElt) -> RngMPolElt {Compte the Hessian of a ternary form.} P3 := Parent(F); require Rank(P3) eq 3: "TernaryHessian(",F,"): Form must be ternary."; n := TotalDegree(F); require forall{m : m in Monomials(F) | TotalDegree(m) eq n}: "TernaryHessian(",F,"): Polynomial must be homogeneous."; dd := func; dF := [dd(F,1,1), dd(F,2,2), dd(F,3,3), dd(F,2,3), dd(F,3,1), dd(F,1,2)]; return dF[1]*dF[2]*dF[3] + 2*dF[4]*dF[5]*dF[6] - dF[1]*dF[4]^2 - dF[2]*dF[5]^2 - dF[3]*dF[6]^2; end intrinsic; P10<[a]> := PolynomialRing(Integers(), 10); S := -( -144*a[1]*a[4]*a[8]*a[10] + 48*a[1]*a[4]*a[9]^2 + 216*a[1]*a[5]*a[7]*a[10] - 24*a[1]*a[5]*a[8]*a[9] - 144*a[1]*a[6]*a[7]*a[9] + 48*a[1]*a[6]*a[8]^2 + 48*a[2]^2*a[8]*a[10] - 16*a[2]^2*a[9]^2 - 144*a[2]*a[3]*a[7]*a[10] + 16*a[2]*a[3]*a[8]*a[9] - 24*a[2]*a[4]*a[5]*a[10] + 16*a[2]*a[4]*a[6]*a[9] + 8*a[2]*a[5]^2*a[9] - 24*a[2]*a[5]*a[6]*a[8] + 48*a[2]*a[6]^2*a[7] + 48*a[3]^2*a[7]*a[9] - 16*a[3]^2*a[8]^2 + 48*a[3]*a[4]^2*a[10] - 24*a[3]*a[4]*a[5]*a[9] + 16*a[3]*a[4]*a[6]*a[8] + 8*a[3]*a[5]^2*a[8] - 24*a[3]*a[5]*a[6]*a[7] - 16*a[4]^2*a[6]^2 + 8*a[4]*a[5]^2*a[6] - a[5]^4); T := 5832*a[1]^2*a[7]^2*a[10]^2 - 3888*a[1]^2*a[7]*a[8]*a[9]*a[10] + 864*a[1]^2*a[7]*a[9]^3 + 864*a[1]^2*a[8]^3*a[10] - 216*a[1]^2*a[8]^2*a[9]^2 - 3888*a[1]*a[2]*a[4]*a[7]*a[10]^2 + 1296*a[1]*a[2]*a[4]*a[8]*a[9]*a[10] - 288*a[1]*a[2]*a[4]*a[9]^3 + 1296*a[1]*a[2]*a[5]*a[7]*a[9]*a[10] - 864*a[1]*a[2]*a[5]*a[8]^2*a[10] + 144*a[1]*a[2]*a[5]*a[8]*a[9]^2 + 1296*a[1]*a[2]*a[6]*a[7]*a[8]*a[10] - 864*a[1]*a[2]*a[6]*a[7]*a[9]^2 + 144*a[1]*a[2]*a[6]*a[8]^2*a[9] + 1296*a[1]*a[3]*a[4]*a[7]*a[9]*a[10] - 864*a[1]*a[3]*a[4]*a[8]^2*a[10] + 144*a[1]*a[3]*a[4]*a[8]*a[9]^2 + 1296*a[1]*a[3]*a[5]*a[7]*a[8]*a[10] - 864*a[1]*a[3]*a[5]*a[7]*a[9]^2 + 144*a[1]*a[3]*a[5]*a[8]^2*a[9] - 3888*a[1]*a[3]*a[6]*a[7]^2*a[10] + 1296*a[1]*a[3]*a[6]*a[7]*a[8]*a[9] - 288*a[1]*a[3]*a[6]*a[8]^3 + 864*a[1]*a[4]^3*a[10]^2 - 864*a[1]*a[4]^2*a[5]*a[9]*a[10] - 864*a[1]*a[4]^2*a[6]*a[8]*a[10] + 576*a[1]*a[4]^2*a[6]*a[9]^2 + 648*a[1]*a[4]*a[5]^2*a[8]*a[10] + 72*a[1]*a[4]*a[5]^2*a[9]^2 + 1296*a[1]*a[4]*a[5]*a[6]*a[7]*a[10] - 720*a[1]*a[4]*a[5]*a[6]*a[8]*a[9] - 864*a[1]*a[4]*a[6]^2*a[7]*a[9] + 576*a[1]*a[4]*a[6]^2*a[8]^2 - 540*a[1]*a[5]^3*a[7]*a[10] - 36*a[1]*a[5]^3*a[8]*a[9] + 648*a[1]*a[5]^2*a[6]*a[7]*a[9] + 72*a[1]*a[5]^2*a[6]*a[8]^2 - 864*a[1]*a[5]*a[6]^2*a[7]*a[8] + 864*a[1]*a[6]^3*a[7]^2 + 864*a[2]^3*a[7]*a[10]^2 - 288*a[2]^3*a[8]*a[9]*a[10] + 64*a[2]^3*a[9]^3 - 864*a[2]^2*a[3]*a[7]*a[9]*a[10] + 576*a[2]^2*a[3]*a[8]^2*a[10] - 96*a[2]^2*a[3]*a[8]*a[9]^2 - 216*a[2]^2*a[4]^2*a[10]^2 + 144*a[2]^2*a[4]*a[5]*a[9]*a[10] + 144*a[2]^2*a[4]*a[6]*a[8]*a[10] - 96*a[2]^2*a[4]*a[6]*a[9]^2 + 72*a[2]^2*a[5]^2*a[8]*a[10] - 48*a[2]^2*a[5]^2*a[9]^2 - 864*a[2]^2*a[5]*a[6]*a[7]*a[10] + 144*a[2]^2*a[5]*a[6]*a[8]*a[9] + 576*a[2]^2*a[6]^2*a[7]*a[9] - 216*a[2]^2*a[6]^2*a[8]^2 - 864*a[2]*a[3]^2*a[7]*a[8]*a[10] + 576*a[2]*a[3]^2*a[7]*a[9]^2 - 96*a[2]*a[3]^2*a[8]^2*a[9] + 144*a[2]*a[3]*a[4]^2*a[9]*a[10] - 720*a[2]*a[3]*a[4]*a[5]*a[8]*a[10] + 144*a[2]*a[3]*a[4]*a[5]*a[9]^2 + 1296*a[2]*a[3]*a[4]*a[6]*a[7]*a[10] - 48*a[2]*a[3]*a[4]*a[6]*a[8]*a[9] + 648*a[2]*a[3]*a[5]^2*a[7]*a[10] - 24*a[2]*a[3]*a[5]^2*a[8]*a[9] - 720*a[2]*a[3]*a[5]*a[6]*a[7]*a[9] + 144*a[2]*a[3]*a[5]*a[6]*a[8]^2 + 144*a[2]*a[3]*a[6]^2*a[7]*a[8] + 144*a[2]*a[4]^2*a[5]*a[6]*a[10] - 96*a[2]*a[4]^2*a[6]^2*a[9] - 36*a[2]*a[4]*a[5]^3*a[10] - 24*a[2]*a[4]*a[5]^2*a[6]*a[9] + 144*a[2]*a[4]*a[5]*a[6]^2*a[8] - 288*a[2]*a[4]*a[6]^3*a[7] + 12*a[2]*a[5]^4*a[9] - 36*a[2]*a[5]^3*a[6]*a[8] + 72*a[2]*a[5]^2*a[6]^2*a[7] + 864*a[3]^3*a[7]^2*a[10] - 288*a[3]^3*a[7]*a[8]*a[9] + 64*a[3]^3*a[8]^3 + 576*a[3]^2*a[4]^2*a[8]*a[10] - 216*a[3]^2*a[4]^2*a[9]^2 - 864*a[3]^2*a[4]*a[5]*a[7]*a[10] + 144*a[3]^2*a[4]*a[5]*a[8]*a[9] + 144*a[3]^2*a[4]*a[6]*a[7]*a[9] - 96*a[3]^2*a[4]*a[6]*a[8]^2 + 72*a[3]^2*a[5]^2*a[7]*a[9] - 48*a[3]^2*a[5]^2*a[8]^2 + 144*a[3]^2*a[5]*a[6]*a[7]*a[8] - 216*a[3]^2*a[6]^2*a[7]^2 - 288*a[3]*a[4]^3*a[6]*a[10] + 72*a[3]*a[4]^2*a[5]^2*a[10] + 144*a[3]*a[4]^2*a[5]*a[6]*a[9] - 96*a[3]*a[4]^2*a[6]^2*a[8] - 36*a[3]*a[4]*a[5]^3*a[9] - 24*a[3]*a[4]*a[5]^2*a[6]*a[8] + 144*a[3]*a[4]*a[5]*a[6]^2*a[7] + 12*a[3]*a[5]^4*a[8] - 36*a[3]*a[5]^3*a[6]*a[7] + 64*a[4]^3*a[6]^3 - 48*a[4]^2*a[5]^2*a[6]^2 + 12*a[4]*a[5]^4*a[6] - a[5]^6; Delta := ExactQuotient(S^3 - T^2, 1728); intrinsic T3Invariants(F::RngMPolElt) -> RngElt, RngElt, RngElt {Compute the invariants S, T and Delta of a ternary cubic.} P3 := Parent(F); require Rank(P3) eq 3: "T3Invariants(",F,"): Form must be ternary."; require forall{m : m in Monomials(F) | TotalDegree(m) eq 3}: "T3Invariants(",F,"): Polynomial must be a cubic form."; mons := [P3.1^i*P3.2^j*P3.3^(3-i-j) : j in [3-i..0 by -1], i in [3..0 by -1]]; cofs := [MonomialCoefficient(F, m) : m in mons]; SF := Evaluate(S, cofs); TF := Evaluate(T, cofs); if Characteristic(P3) notin {2, 3} then DF := ExactQuotient(SF^3 - TF^2, Parent(SF)!1728); else DF := Evaluate(Delta, cofs); end if; return SF, TF, DF; end intrinsic; ///////////////////////////////////////////////////////////////////// // Minimisation of cubics ///////////////////////////////////////////////////////////////////// intrinsic Minimise(C::RngMPolElt) -> RngMPolElt, AlgMatElt {Minimise the integral ternary cubic form C. Returns the minimised form and the transformation.} P3 := Parent(C); require Rank(P3) eq 3: "Form must be ternary."; require forall{m : m in Monomials(C) | TotalDegree(m) eq 3}: "Polynomial must be a cubic form."; require BaseRing(P3) cmpeq Rationals(): "Form must live in polynomial ring over the rationals."; require forall{c : c in Coefficients(C) | Denominator(c) eq 1}: "Form must have integral coefficients."; vprintf Minimise, 1: "Minimise C =\n%o\n", C; trans := IdentityMatrix(Rationals(), 3); S, T, D := T3Invariants(C); g := GCD([Integers() | S, T, D]); vprintf Minimise, 2: " Factoring GCD of invariants ...\n"; plist := [p : p in PrimeDivisors(g) | (S eq 0 or Valuation(S, p) ge 4) and (T eq 0 or Valuation(T, p) ge 6) and (D eq 0 or Valuation(D, p) ge 12)]; vprintf Minimise, 2: " primes where model is nonminimal:\n %o\n", plist; for p in plist do C, tr := pMinimise(C, p); trans := trans*tr; end for; vprintf Minimise, 1: "Minimise -- done.\n"; return C, trans; end intrinsic; intrinsic pMinimise(C::RngMPolElt, p::RngIntElt) -> RngMPolElt, AlgMatElt {Minimise the integral ternary cubic Form C at the prime p. Returns the minimised form and the transformation.} P3 := Parent(C); require Rank(P3) eq 3: "Form must be ternary."; require forall{m : m in Monomials(C) | TotalDegree(m) eq 3}: "Polynomial must be a cubic form."; require BaseRing(P3) cmpeq Rationals(): "Form must live in polynomial ring over the rationals."; require forall{c : c in Coefficients(C) | Denominator(c) eq 1}: "Form must have integral coefficients."; vprintf Minimise, 1: "pMinimise at p = %o ...\n", p; C0 := C; // save for later check trans := IdentityMatrix(Rationals(), 3); indices := [ : j in [0..3-i], i in [0..3]]; mons := [P3.1^(a[1])*P3.2^(a[2])*P3.3^(a[3]) : a in indices]; MC := MonomialCoefficient; F := GF(p); PF := PolynomialRing(F, 3); S, T, D := T3Invariants(C); steps := Min([Floor(Valuation(pair[1],p)/pair[2]) : pair in [, , ] | pair[1] ne 0]); if p eq 3 then // check Kraus's conditions Tred := Integers()!(T/3^(6*steps)); if Tred mod 27 in {9, 18} then steps -:= 1; end if; elif p eq 2 then // check Kraus's conditions Tred := Integers()!(T/2^(6*steps)); if Tred mod 4 ne 3 and (Valuation(S, 2) lt 4*steps+4 or Tred mod 32 notin {0, 8}) then steps -:= 1; end if; end if; // now steps is the number of minimisation steps to make vprintf Minimise, 2: " %o minimisation steps to make.\n", steps; while steps gt 0 do vprintf Minimise, 3: " remaining steps: %o\n", steps; CF := PF!C; if CF eq 0 then C := 1/p*C; vprintf Minimise, 2: " C mod %o is zero. Dividing C by %o.\n", p, p; steps -:= 1; else fact := Factorisation(CF); if forall(f){f : f in fact | f[2] eq 1} then vprintf Minimise, 3: " C mod %o is squarefree.\n", p; I := ideal; assert Dimension(I) eq 1; // isolated singularities // we know that there is one singular point; move it to (0:0:1) V1 := Variety(ideal); if IsEmpty(V1) then V2 := Variety(ideal); if IsEmpty(V2) then V3 := Variety(ideal); assert not IsEmpty(V3); sing := V3[1]; // = (1:0:0) C := Evaluate(C, [P3.3, P3.2, -P3.1]); trans := trans*Matrix(Rationals(), [[0,0,1], [0,1,0], [-1,0,0]]); else // sing = (x:1:0) sing := V2[1]; C := Evaluate(C, [P3.1 + (Integers()!sing[1])*P3.3, P3.3, -P3.2]); trans := trans*Matrix(Rationals(), [[1,0,Integers()!sing[1]], [0,0,1], [0,-1,0]]); end if; else sing := V1[1]; C := Evaluate(C, [P3.1 + (Integers()!sing[1])*P3.3, P3.2 + (Integers()!sing[2])*P3.3, P3.3]); trans := trans*Matrix(Rationals(), [[1,0,Integers()!sing[1]], [0,1,Integers()!sing[2]], [0,0,1]]); end if; CF := PF!C; // look at degree 2 part in X,Y deg2 := Evaluate(CF, [PF.1, PF.2, 1]) - Evaluate(CF, [PF.1, PF.2, 0]); if deg2 eq 0 then // CF only depends on X,Y, so describes three lines through one point // since CF is squarefree, the lines are distinct // if minimisation is possible, then by z --> z/p vprintf Minimise, 3: " C mod %o is union of three distinct lines.\n", p; if forall{i : i in [1..#indices] | Valuation(MC(C, mons[i]), p) ge indices[i,3]} then C := Evaluate(C, [P3.1, P3.2, 1/p*P3.3]); trans := trans*Matrix(Rationals(), [[1,0,0], [0,1,0], [0,0,1/p]]); steps -:= 1; else break; end if; else // have nonzero degree 2 part // degree 2 part must be a square vprintf Minimise, 3: " C mod %o is cuspidal cubic or line touching conic.\n", p; a := MC(deg2, PF.1^2); b := MC(deg2, PF.1*PF.2); c := MC(deg2, PF.2^2); assert 4*a*c eq b^2; // transform to Y^2 if c ne 0 then if p eq 2 then s := a/c; else s := b/(2*c); end if; C := Evaluate(C, [P3.1, P3.2 - (Integers()!s)*P3.1, P3.3]); trans := trans*Matrix(Rationals(), [[1,0,0], [-Integers()!s,1,0], [0,0,1]]); else // c eq 0 C := Evaluate(C, [P3.2, -P3.1, P3.3]); trans := trans*Matrix(Rationals(), [[0,1,0], [-1,0,0], [0,0,1]]); end if; // c ne 0 CF := PF!C; // now CF = c*Y^2*Z + F(X,Y) must be either a cuspidal cubic // or a line touching a conic; in both cases we must have // p^2 | a003, a102 assert Valuation(MC(C, P3.3^3), p) ge 2 and Valuation(MC(C, P3.1*P3.3^2), p) ge 2; vprintf Minimise, 3: " transform to another model of same level.\n"; C := p*Evaluate(C, [P3.1, P3.2, 1/p*P3.3]); trans := trans*Matrix(Rationals(), [[1,0,0], [0,1,0], [0,0,1/p]]); // this leads to CF being a product of linear factors end if; // deg2 eq 0 else // a multiple line // f = // move multiple line to Y line := f[1]; a,b,c := Explode([MC(line, PF.i) : i in [1..3]]); if b ne 0 then C := Evaluate(C, [P3.1, P3.2 - (Integers()!(a/b))*P3.1 - (Integers()!(c/b))*P3.3, P3.3]); trans := trans*Matrix(Rationals(), [[1,0,0], [-Integers()!(a/b),1,-Integers()!(c/b)], [0,0,1]]); elif c ne 0 then // line = aX + cZ C := Evaluate(C, [P3.1, -P3.3, P3.2 - (Integers()!(a/c))*P3.1]); trans := trans*Matrix(Rationals(), [[1,0,0], [0,0,-1], [-Integers()!(a/c),1,0]]); else // line = aX C := Evaluate(C, [P3.2, -P3.1, P3.3]); trans := trans*Matrix(Rationals(), [[0,1,0], [-1,0,0], [0,0,1]]); end if; CF := PF!C; if f[2] eq 2 then // a double line and a simple one // move to Y^2 Z vprintf Minimise, 3: " C mod %o is the union of a double line and a simple line.\n", p; line := ExactQuotient(CF, PF.2^2); a,b,c := Explode([MC(line, PF.i) : i in [1..3]]); if c ne 0 then C := Evaluate(C, [P3.1, P3.2, P3.3 - (Integers()!(a/c))*P3.1 - (Integers()!(b/c))*P3.2]); trans := trans*Matrix(Rationals(), [[1,0,0], [0,1,0], [-Integers()!(a/c),-Integers()!(b/c),1]]); else assert a ne 0; C := Evaluate(C, [P3.3, P3.2, -P3.1]); trans := trans*Matrix(Rationals(), [[0,0,1], [0,1,0], [-1,0,0]]); end if; CF1 := PF!(1/p*&+[MC(C, P3.1^i*P3.3^(3-i))*P3.1^i*P3.3^(3-i) : i in [0..3]]); if MC(CF1, PF.1^3) eq 0 then if forall{i : i in [1..#indices] | Valuation(MC(C, mons[i]), p) ge indices[i,1]} then C := Evaluate(C, [1/p*P3.1, P3.2, P3.3]); trans := trans*Matrix(Rationals(), [[1/p,0,0], [0,1,0], [0,0,1]]); steps -:= 1; elif forall{i : i in [1..#indices] | Valuation(MC(C, mons[i]), p) ge indices[i,1]+indices[i,3]-1} then C := p*Evaluate(C, [1/p*P3.1, P3.2, 1/p*P3.3]); trans := trans*Matrix(Rationals(), [[1/p,0,0], [0,1,0], [0,0,1/p]]); steps -:= 1; else break; end if; else // MC(CF1, PF.1^3) ne 0 // CF1 must be a cube fact := Factorisation(CF1); if fact[1,2] ne 3 then break; end if; r := -MC(fact[1,1], PF.3)/MC(fact[1,1], PF.1); // triple root C := Evaluate(C, [P3.1 + (Integers()!r)*P3.3, P3.2, P3.3]); trans := trans*Matrix(Rationals(), [[1,0,Integers()!r], [0,1,0], [0,0,1]]); if forall{i : i in [1..#indices] | Valuation(MC(C, mons[i]), p) ge 1-indices[i,2]} then vprintf Minimise, 3: " transform to another model of same level.\n"; C := 1/p*Evaluate(C, [P3.1, p*P3.2, P3.3]); trans := trans*Matrix(Rationals(), [[1,0,0], [0,p,0], [0,0,1]]); // makes CF a triple line else break; end if; end if; // MC(CF1, PF.1^3) eq 0 else // a triple line // CF = Y^3 vprintf Minimise, 3: " C mod %o is a triple line.\n", p; CF1 := PF!(1/p*&+[MC(C, P3.1^i*P3.3^(3-i))*P3.1^i*P3.3^(3-i) : i in [0..3]]); // this must have a multiple root if CF1 eq 0 then C := 1/p^2*Evaluate(C, [P3.1, p*P3.2, P3.3]); trans := trans*Matrix(Rationals(), [[1,0,0], [0,p,0], [0,0,1]]); steps -:= 1; else fact := Factorisation(CF1); if not exists(f){f : f in fact | f[2] ge 2} then break; end if; lin := f[1]; a := MC(lin, PF.1); b := MC(lin, PF.3); if a ne 0 then C := Evaluate(C, [P3.1 - (Integers()!(b/a))*P3.3, P3.2, P3.3]); trans := trans*Matrix(Rationals(), [[1,0,-Integers()!(b/a)], [0,1,0], [0,0,1]]); else C := Evaluate(C, [-P3.3, P3.2, P3.1]); trans := trans*Matrix(Rationals(), [[0,0,-1], [0,1,0], [1,0,0]]); end if; if forall{i : i in [1..#indices] | Valuation(MC(C, mons[i]), p) ge indices[i,3]} then C := Evaluate(C, [P3.1, P3.2, 1/p*P3.3]); trans := trans*Matrix(Rationals(), [[1,0,0], [0,1,0], [0,0,1/p]]); steps -:= 1; else break; end if; end if; end if; // if f[2] eq 2 end if; // forall(f){f : f in fact | f[2] eq 1} end if; assert IsDivisibleBy(Evaluate(C0, [&+[trans[i,j]*P3.j : j in [1..3]] : i in [1..3]]), C); end while; if steps gt 0 then vprintf Minimise, 1: "pMinimise at p = %o: not completely minimised.\n", p; end if; vprintf Minimise, 1: "pMinimise at p = %o -- done.\n", p; l := LCM([Denominator(c) : c in Eltseq(trans)]); trans *:= l; g := GCD([Numerator(c) : c in Eltseq(trans)]); trans *:= 1/g; return C, trans; end intrinsic; // Convert (E, pt) --> (C, pt') where E is an elliptic curve/Q with // a rational point pt and C is the corresponding plane cubic, with a // rational point pt' that maps to pt on E. intrinsic CubicFromPoint(E::CrvEll, P::PtEll) -> RngMPolElt, SeqEnum {Given an elliptic curve and a point on it, find the corresponding 3-covering as a ternary cubic, together with a point on it that maps to the given point on E.} // Bring (E, P) into the form // E : y(y-b) = (x-a)(x^2+c), P = (a,0) R := CoefficientRing(E); require Characteristic(R) ne 2: "This does not work in characteristic 2."; a1,a2,a3,a4,a6 := Explode(aInvariants(E)); require P[3] ne 0: "Point may not be at infinity."; px := P[1]/P[3]; py := P[2]/P[3]; P := PolynomialRing(R, 2); pol := Y^2 + a1*X*Y + a3*Y - (X^3 + a2*X^2 + a4*X + a6); // shift point to (0,0) pol := Evaluate(pol, [X+px, Y+py]); assert Evaluate(pol, [0,0]) eq 0; // get rid of the xy term pol := Evaluate(pol, [X, Y-1/(R!2)*MonomialCoefficient(pol,X*Y)*X]); assert MonomialCoefficient(pol, X*Y) eq 0 and Evaluate(pol, [0,0]) eq 0; // shift x to get RHS into form (x-a)(x^2+c) (from x(x^2+a2*x+a4)) a := -1/(R!2)*MonomialCoefficient(pol, X^2); pol := Evaluate(pol, [X-a, Y]); assert MonomialCoefficient(ExactQuotient(Coefficient(pol,Y,0),(X-a)),X) eq 0; assert Evaluate(pol, [a,0]) eq 0; b := -MonomialCoefficient(pol,Y); c := -MonomialCoefficient(pol,X); P3 := PolynomialRing(R, 3); // cubic and point C := y*z^2 + x^2*z + a*x^2*y + b*x*y^2 + c*y^3; ps := [R | 0,0,1]; // if R = rationals, then minimise and reduce if R cmpeq Rationals() then l := LCM([Denominator(a), Denominator(b), Denominator(c)]); C *:= l; C, tr := Minimise(C); ps := Eltseq(tr^-1*Matrix(R, [[c] : c in ps])); C, tr := Reduce(C); ps := Eltseq(ChangeRing(tr^-1,R)*Matrix(R, [[c] : c in ps])); l := LCM([Denominator(c) : c in ps]); ps := [l*c : c in ps]; g := GCD([Numerator(c) : c in ps]); ps := [c/g : c in ps]; assert Evaluate(C, ps) eq 0; end if; return C, ps; end intrinsic; /* // Idea: Given that we know by how much we have to minimise the given // cubic (it is N = (disc(cubic)/minimal_disc(Jac(cubic)))^(1/12)), // we just enter into the computation (excepting p = 2 and 3). If // something goes wrong, then because we have an element mod N that // is neither zero nor invertible and thus will split N. In this case, // split N and repeat with the factors. forward extendBD; function split(a,b) // printf "split(%o,%o)\n", a, b; g := GCD(a,b); if g eq 1 then return [c : c in [a,b] | c ne 1]; end if; a := ExactQuotient(a,g); b := ExactQuotient(b,g); return extendBD([a,b], g); end function; function extendBD(l, n) // printf "extendBD(%o,%o)\n", l, n; gl := [GCD(a, n) : a in l]; l1 := [l[i] : i in [1..#l] | gl[i] eq 1]; l2 := [l[i] : i in [1..#l] | gl[i] ne 1]; if IsEmpty(l2) then return [c : c in l cat [n] | c ne 1]; end if; s := split(l2[1], n); res := l2[2..#l2]; for a in s do res := extendBD(res, a); end for; return l1 cat res; end function; function BD(l) res := []; for a in l do res := extendBD(res, a); end for; return res; end function; function BreakUpWith(n, l) // use the numbers in the sequence l to split n into factors bl := BD([n] cat [GCD(n,a) : a in l]); function val(m, p) e := 0; while IsDivisibleBy(m, p) do e +:= 1; m := m div p; end while; return m, e; end function; res := []; for b in bl do n, e := val(n, b); Append(~res, ); end for; return res; end function; intrinsic Minimise(C::RngMPolElt, l::SeqEnum : Check := true) -> RngMPolElt, Mtrx {Minimise C according to the factorisation sequence l (the factors need not to be prime).} vprintf Minimise, 1: "Minimise(l = %o)\n", l; tr := IdentityMatrix(Rationals(), 3); for ll in l do for i := 1 to ll[2] do C, tr1 := Minimise(C, ll[1] : Check := Check); tr := tr*tr1; end for; end for; return C, tr; end intrinsic; //-------------------------------------------------------------- function val(m, p) if m eq 0 then return Infinity(); end if; m := Integers()!m; p := Integers()!p; e := 0; while IsDivisibleBy(m ,p) do e +:= 1; m := m div p; end while; return e; end function; function valmax(m, p, max) m := Integers()!m; p := Integers()!p; e := 0; while IsDivisibleBy(m ,p) do e +:= 1; m := m div p; if e ge max then return e; end if; end while; return e; end function; intrinsic MinimiseQuick(C::RngMPolElt) -> RngMPolElt {Do all minimisation that is possible by scalings onnly.} P3 := Parent(C); require Rank(P3) eq 3: "Form must be ternary."; require forall{m : m in Monomials(C) | TotalDegree(m) eq 3}: "Polynomial must be a cubic form."; require BaseRing(P3) cmpeq Rationals(): "Form must live in polynomial ring over the rationals."; require forall{c : c in Coefficients(C) | Denominator(c) eq 1}: "Form must have integral coefficients."; C := 1/GCD(ChangeUniverse(Coefficients(C), Integers()))*C; S, T, D := T3Invariants(C); S := Integers()!S; T := Integers()!T; D := Integers()!D; list := [l[1] : l in BreakUpWith(GCD([S, T, D]), ChangeUniverse(Coefficients(C), Integers()))]; Sort(~list); Reverse(~list); mons := MonomialsOfDegree(P3, 3); for n in list do repeat vals := [valmax(MonomialCoefficient(C, m), n, 10) : m in mons]; inds := [ : m in mons]; // find best possible scaling test := func; test1 := function(sy) sz := 0; t0 := test(sy,0); tp := test(sy,1); tm := test(sy,-1); if tp gt t0 then sz := 1; repeat sz +:= 1; t0 := tp; tp := test(sy,sz); until tp le t0; return sz-1, t0; elif tm gt t0 then sz := -1; repeat sz -:= 1; t0 := tm; tm := test(sy,sz); until tm le t0; return sz+1, t0; else return 0, t0; end if; end function; sy := 0; sz,t0 := test1(0); szp,tp := test1(1); szm,tm := test1(-1); if tp gt t0 then sy := 1; repeat sy +:= 1; t0 := tp; sz := szp; szp; tp := test1(sy); until tp le t0; sy -:= 1; elif tm gt t0 then sy := -1; repeat sy -:= 1; t0 := tm; sz := szm; szm, tm := test1(sy); until tm le t0; sy +:= 1; end if; C := n^-(t0+sy+sz)*Evaluate(C, [P3.1, n^sy*P3.2, n^sz*P3.3]); assert forall{c : c in Coefficients(C) | Denominator(c) eq 1}; until t0 eq 0; end for; return C; end intrinsic; //-------------------------------------------------------------- function Mini300(C, N) // here C mod N is const*x^3 printf "Not yet implemented (Mini300) -- sorry.\n"; return C; end function; function Mini210(C, N) // here C mod N is const*x^2*y vprintf Minimise, 2: "Mini210(N = %o)\n", N; return C; end function; function Mini2var(C, N) // here C mod N does not involve z F := Integers(N); PF := PolynomialRing(F,3); Cred := PF!C; vprintf Minimise, 2: "Mini2Var(Cred= %o, N = %o)\n", Cred, N; lc := MonomialCoefficient(Cred, PF.1^3); if lc eq 0 then // C = y (a x^2 + b x y + c y^2) mod N lc := MonomialCoefficient(Cred, PF.1^2*PF.2); if lc eq 0 then // C = y^2 (a x + b y) mod N lc := MonomialCoefficient(Cred, PF.1*PF.2^2); if lc eq 0 then // C = c y^3 mod N lc := MonomialCoefficient(Cred, PF.2^3); assert IsUnit(lc); tr1 := Matrix([[0,-1,0], [1,0,0], [0,0,1]]); return Mini300(PolMat(C, tr1), N); elif IsUnit(lc) then // x <- -y, y <- x - b/a y b := MonomialCoefficient(Cred, PF.2^3); tr1 := Matrix([[0,-1,0], [1,-Integers()!(b/lc),0], [0,0,1]]); return Mini210(PolMat(C, tr1), N); else return Minimise(C, BreakUpWith(N, [Integers() | lc]) : Check := false); end if; elif IsUnit(lc) then // check discriminant b := MonomialCoefficient(Cred, PF.1*PF.2^2); c := MonomialCoefficient(Cred, PF.2^3); disc := b^2 - 4*lc*c; if disc eq 0 then // C = y (...)^2 mod N // find linear factor whose square we know if IsEven(N) and N ne 2 then return Minimise(C, BreakUpWith(N, [Integers() | 2]) : Check := false); elif N eq 2 then // b = 0 mod 2, factor is lc*x + c*y = x + c*y tr1 := Matrix([[1,-Integers()!(c/lc),0], [0,1,0], [0,0,1]]); return Mini210(PolMat(C, tr1), N); else // N is odd -> get x + b/(2*lc)*y tr1 := Matrix([[1,-Integers()!(b/(2*lc)),0], [0,1,0], [0,0,1]]); return Mini210(PolMat(C, tr1), N); end if; elif IsUnit(disc) then // three distinct roots Cnew := PolMat(C, DiagonalMatrix([1,1,1/N])); if forall{c : c in Coefficients(Cnew) | Denominator(c) eq 1} then return Cnew; else printf "Minimisation at N = %o not possible.\n", N; return C; end if; else return Minimise(C, BreakUpWith(N, [Integers() | disc]) : Check := false); end if; else return Minimise(C, BreakUpWith(N, [Integers() | lc]) : Check := false); end if; elif IsUnit(lc) then P1 := PolynomialRing(F); C1 := Evaluate(Cred, [P1.1,1,0]); a := Coefficient(C1, 2); b := Coefficient(C1, 1); c := Coefficient(C1, 0); disc := Discriminant(C1); if disc eq 0 then // multiple root DC1 := Derivative(C1); if DC1 eq 0 then // must have N eq 3 assert N eq 3; tr1 := Matrix([[1,-Integers()!(c/lc),0], [0,1,0], [0,0,1]]); return Mini300(PolMat(C, tr1), N); end if; disc1 := Discriminant(DC1); if disc1 eq 0 then // triple root if N eq 3 then // -> (x + c/lc*y)^3 tr1 := Matrix([[1,-Integers()!(c/lc),0], [0,1,0], [0,0,1]]); return Mini300(PolMat(C, tr1), N); elif IsDivisibleBy(N, 3) then return Minimise(C, BreakUpWith(N, [3]) : Check := false); else // N prime to 3 -> (x + a/(3*lc)*y)^3 tr1 := Matrix([[1,-Integers()!(a/(3*lc)),0], [0,1,0], [0,0,1]]); return Mini300(PolMat(C, tr1), N); end if; elif IsUnit(disc1) then if N eq 3 then // DC1 = 2*a*u + b assert IsUnit(a); dr := -b/(2*a); // double root sr := -(a/lc + 2*dr); // single root tr1 := Matrix([[1,Integers()!dr,0], [0,1,0], [0,0,1]]); sr := sr-dr; assert IsUnit(sr); tr1 := tr1*Matrix([[1,0,0], [Integers()!(1/sr),1,0], [0,0,1]]); return Mini210(PolMat(C, tr1), N); elif IsDivisibleBy(N, 3) then return Minimise(C, BreakUpWith(N, [3]) : Check := false); else DC1 := 1/LeadingCoefficient(DC1)*DC1; lin := C1 - lc*P1.1*DC1; lin := lin - Coefficient(lin, 2)*DC1; ll := Coefficient(lin ,1); assert IsUnit(ll); dr := -Coefficient(lin,0)/ll; // double root sr := -(a/lc + 2*dr); // single root tr1 := Matrix([[1,Integers()!dr,0], [0,1,0], [0,0,1]]); sr := sr-dr; assert IsUnit(sr); tr1 := tr1*Matrix([[1,0,0], [Integers()!(1/sr),1,0], [0,0,1]]); return Mini210(PolMat(C, tr1), N); end if; else return Minimise(C, BreakUpWith(N, [Integers() | disc1]) : Check := false); end if; elif IsUnit(disc) then // three distinct roots Cnew := PolMat(C, DiagonalMatrix([1,1,1/N])); if forall{c : c in Coefficients(Cnew) | Denominator(c) eq 1} then return Cnew; else printf "Minimisation at N = %o not possible.\n", N; return C; end if; else return Minimise(C, BreakUpWith(N, [Integers() | disc]) : Check := false); end if; else return Minimise(C, BreakUpWith(N, [Integers() | lc]) : Check := false); end if; return C; end function; function MinLine(C, N, line) vprintf Minimise, 2: " MinLine(N = %o, line = %o)\n", N, line; F := Integers(N); PF := PolynomialRing(F,3); if line[1] eq 1 then line[2] := -line[2]; line[3] :=-line[3]; tr1 := Matrix(Rationals(), [ChangeUniverse(line, Integers()), [0,1,0], [0,0,1]]); elif line[2] eq 1 then line[1] := -line[1]; line[3] :=-line[3]; tr1 := Matrix(Rationals(), [[1,0,0], ChangeUniverse(line, Integers()), [0,0,1]]); tr1 := tr1*Matrix(Rationals(), [[0,1,0], [-1,0,0], [0,0,1]]); elif line[3] eq 1 then line[1] := -line[1]; line[2] :=-line[2]; tr1 := Matrix(Rationals(), [[1,0,0], [0,1,0], ChangeUniverse(line, Integers())]); tr1 := tr1*Matrix(Rationals(), [[0,0,1], [0,1,0], [-1,0,0]]); else error "line equation must have one coefficient = 1."; end if; C := PolMat(C, tr1); Cred := PF!C; assert Setseq(Monomials(Cred)) subset {PF.1^2*PF.i : i in [1..3]}; fact2 := &+[PF.i*MonomialCoefficient(Cred, PF.1^2*PF.i) : i in [1..3]]; a,b,c := Explode([MonomialCoefficient(fact2, PF.i) : i in [1..3]]); if b eq 0 and c eq 0 then return Mini300(C, N); elif IsUnit(b) then tr1 := Matrix(Rationals(), [[1,0,0], [Integers() | -a/b, 1, -c/b], [0,0,1]]); return Mini210(PolMat(C, tr1), N); elif IsUnit(c) then tr1 := Matrix(Rationals(), [[1,0,0], [0,1,0], [Integers() | -a/c, -b/c, 1]]); return Mini210(PolMat(C, tr1), N); else return Minimise(C, BreakUpWith(N, [Integers() | b,c]) : Check := false); end if; end function; function MinSing1(C, N) // here C = c x^2 z + F(x,y) mod N, c invertible vprintf Minimise, 2: " MinSing1(N = %o)\n", N; F := Integers(N); PF := PolynomialRing(F,3); Cred := PF!C; g := GCD([Integers() | MonomialCoefficient(Cred, PF.1^i*PF.2^(3-i)) : i in [0..3]]); if g eq 1 then P3 := Parent(C); Cnew := N*Evaluate(C, [P3.1, P3.2, 1/N*P3.3]); if forall{c : c in Coefficients(Cnew) | Denominator(c) eq 1} then tr1 := Matrix([[1,0,0], [0,0,-1], [0,1,0]]); return Mini2var(PolMat(Cnew, tr1), N); else printf "Minimisation not possible.\n"; return C; end if; elif F!g eq 0 then tr1 := Matrix([[1,0,0], [0,0,-1], [0,1,0]]); return Mini210(PolMat(C, tr1), N); else return Minimise(C, BreakUpWith(N, [g]) : Check := false); end if; end function; function MinSing(C, N, sing) // here sing are coordinates of a singular point mod N vprintf Minimise, 2: " MinSing(N = %o, sing = %o)\n", N, sing; F := Integers(N); PF := PolynomialRing(F,3); // move sing to (0 : 0 : 1) if sing[1] eq 1 then tr1 := Matrix([[1,0,0], [Integers() | sing[2],1,0], [Integers() | sing[3],0,1]]); tr1 := tr1*Matrix([[0,0,-1], [0,1,0], [1,0,0]]); elif sing[2] eq 1 then tr1 := Matrix([[Integers() | 1,sing[1],0], [0,1,0], [Integers() | 0,sing[3],1]]); tr1 := tr1*Matrix([[1,0,0], [0,0,-1], [0,1,0]]); elif sing[3] eq 1 then tr1 := Matrix([[Integers() | 1,0,sing[1]], [Integers() | 0,1,sing[2]], [0,0,1]]); else error "One of the coordinates of sing must be 1."; end if; C := PolMat(C, tr1); Cred := PF!C; assert forall{m : m in {PF.3^3, PF.1*PF.3^2, PF.2*PF.3^2} | MonomialCoefficient(Cred, m) eq 0}; // now must have that degree-2 part in (x,y) is a square a, b, c := Explode([MonomialCoefficient(Cred, PF.1^(2-i)*PF.2^i*PF.3) : i in [0..2]]); if b^2 ne 4*a*c then printf "Minimisation not possible.\n"; return C; end if; if a eq 0 then if b ne 0 then return Minimise(C, BreakUpWith(N, [Integers() | b]) : Check := false); elif c eq 0 then // no F_2(x,y)*z return Mini2var(C, N); elif IsUnit(c) then tr1 := Matrix([[0,-1,0], [1,0,0], [0,0,1]]); return MinSing1(PolMat(C, tr1), N); else return Minimise(C, BreakUpWith(N, [Integers() | c]) : Check := false); end if; elif IsUnit(a) then if N eq 2 then tr1 := Matrix([[1,-Integers()!(c/a),0], [0,1,0], [0,0,1]]); return MinSing1(PolMat(C, tr1), N); elif IsDivisibleBy(N, 2) then return Minimise(C, BreakUpWith(N, [2]) : Check := false); else // N is odd tr1 := Matrix([[1,-Integers()!(b/(2*a)),0], [0,1,0], [0,0,1]]); return MinSing1(PolMat(C, tr1), N); end if; else return Minimise(C, BreakUpWith(N, [Integers() | a]) : Check := false); end if; return C; end function; intrinsic Minimise(C::RngMPolElt, N::RngIntElt : Check := true) -> RngMPolElt, Mtrx {Minimise C by N.} vprintf Minimise, 1: "Minimise(N = %o),\n Curve = %o\n", N, C; N := Abs(N); P3 := Parent(C); require Rank(P3) eq 3: "Form must be ternary."; require forall{m : m in Monomials(C) | TotalDegree(m) eq 3}: "Polynomial must be a cubic form."; require BaseRing(P3) cmpeq Rationals(): "Form must live in polynomial ring over the rationals."; require forall{c : c in Coefficients(C) | Denominator(c) eq 1}: "Form must have integral coefficients."; C := 1/GCD(ChangeUniverse(Coefficients(C), Integers()))*C; if Check then S,T := T3Invariants(C); S := Integers()!S; T := Integers()!T; require IsDivisibleBy(S, N^4) and IsDivisibleBy(T, N^6): "N^4, N^6 do not divide the invariants."; end if; Corig := C; if N eq 1 then return C; end if; F := Integers(N); PF := PolynomialRing(F, 3); Cred := PF!C; derivs := [Derivative(Cred, i) : i in [1..3]]; if derivs[1] eq 0 then // no x return Mini2var(C, N); elif derivs[2] eq 0 then // no y tr1 := Matrix(Rationals(), [[0,1,0], [-1,0,0], [0,0,1]]); return Mini2var(PolMat(C, tr1), N); elif derivs[3] eq 0 then // no z tr1 := Matrix(Rationals(), [[0,0,1], [0,1,0], [-1,0,0]]); return Mini2var(PolMat(C, tr1), N); end if; gl := [GCD(N, GCD(ChangeUniverse(Coefficients(Fd), Integers()))) : Fd in derivs]; if exists{g : g in gl | g ne 1} then return Minimise(Corig, BreakUpWith(N, gl) : Check := false); end if; // here, all the derivatives are of content prime to N mat := Matrix([[MonomialCoefficient(Fd, PF.1*PF.j) : j in [1..3]] : Fd in derivs]); det := Determinant(mat); if det eq 0 then // distinguish further depending on rank of mat minors := Minors(mat, 2); if forall{m : m in minors | m eq 0} then // rank is at most 1 if mat eq 0 then // rank is 0 // derivatives do not involve x -- singularities along a line // derivatives must have common factor lcs := [MonomialCoefficient(Fd, PF.2^2) : Fd in derivs]; _, pos := Position([IsUnit(l) : l in lcs], true); if forall{c : c in lcs | c eq 0} then // z is common factor, line is y = 0 return MinLine(C, N, [F | 0, 1, 0]); elif pos ne 0 then q1 := derivs[pos]; qs := [derivs[i] : i in [1..3] | i ne pos]; q2 := qs[1] - MonomialCoefficient(qs[1], PF.2^2)*q1; q2 := &+[MonomialCoefficient(q2, PF.2*m)*m : m in [PF.2,PF.3]]; q3 := qs[2] - MonomialCoefficient(qs[2], PF.2^2)*q1; q3 := &+[MonomialCoefficient(q3, PF.2*m)*m : m in [PF.2,PF.3]]; lq2 := MonomialCoefficient(q2, PF.2); lq3 := MonomialCoefficient(q3, PF.2); assert lq2 ne 0 and lq3 ne 0; if IsUnit(lq2) then return MinLine(C, N, [0, 1, MonomialCoefficient(q2, PF.3)/lq2]); elif IsUnit(lq3) then return MinLine(C, N, [0, 1, MonomialCoefficient(q3, PF.3)/lq3]); else return Minimise(Corig, BreakUpWith(N, [Integers() | lq2, lq3]) : Check := false); end if; else return Minimise(Corig, BreakUpWith(N, ChangeUniverse(lcs, Integers())) : Check := false); end if; // all lcs = 0 elif exists{c : c in Eltseq(mat) | IsUnit(c)} then // rank is 1 _, pos := Position([IsUnit(c) : c in Eltseq(mat)], true); ii := 1 + ((pos-1) div 3); jj := 1 + ((pos-1) mod 3); p1 := 1/mat[ii,jj]*derivs[ii]; ps := [derivs[i] - MonomialCoefficient(derivs[i], PF.1*PF.jj)*p1 : i in [1..3] | i ne ii]; assert forall{p : p in ps | Seqset(Monomials(p)) subset {PF.2^i*PF.3^(2-i) : i in [0..2]}}; // two possibilities: // - either ps[1] = ps[2] = 0, and p1 is a square // - or ps[1] and ps[2] have a common root if ps[1] eq 0 and ps[2] eq 0 then if IsEven(N) and N ne 2 then return Minimise(Corig, BreakUpWith(N, [2]) : Check := false); end if; if N eq 2 then assert forall{m : m in [PF.1*PF.2, PF.2*PF.3, PF.3*PF.1] | MonomialCoefficient(p1, m) eq 0}; return MinLine(C, N, [MonomialCoefficient(p1, PF.i^2) : i in [1..3]]); else // N is odd assert jj eq 1; // so p1 = x^2 + ... assert Coefficient(p1, PF.1, 1)^2 eq 4*Coefficient(p1, PF.1, 0); return MinLine(C, N, [F | 1, MonomialCoefficient(p1, PF.1*PF.2)/2, MonomialCoefficient(p1, PF.1*PF.3)/2]); end if; end if; q2 := ps[1]; q3 := ps[2]; lq2 := MonomialCoefficient(q2, PF.2); lq3 := MonomialCoefficient(q3, PF.2); if lq2 eq 0 and lq3 eq 0 then // z is common factor p1 := Evaluate(p1, [PF.1, PF.2, 0]); if p1 eq 0 then return MinLine(C, N, [F | 0, 0, 1]); end if; g := GCD(N, GCD(ChangeUniverse(Coefficients(p1), Integers()))); if g ne 1 then return Minimise(Corig, BreakUpWith(N, [g]) : Check := false); else // p1 must be a square lp1 := MonomialCoefficient(p1, PF.1^2); if IsEven(N) and N ne 2 then return Minimise(Corig, BreakUpWith(N, [2]) : Check := false); elif IsUnit(lp1) then if N eq 2 then assert MonomialCoefficient(p1, PF.1*PF.2) eq 0; return MinSing(C, N, [1, MonomialCoefficient(p1, PF.2^2), 0]); else a, b, c := Explode([MonomialCoefficient(p1, PF.1^(2-i)*PF.2^i) : i in [0..2]]); assert b^2 eq 4*a*c; return MinSing(C, N, [1, -b/(2*a), 0]); end if; elif lp1 eq 0 then assert MonomialCoefficient(p1, PF.1*PF.2) eq 0; return MinSing(C, N, [F | 1, 0, 0]); else return Minimise(Corig, BreakUpWith(N, [Integers()!lp1]) : Check := false); end if; end if; elif IsUnit(lq2) then // ??? printf "Not yet implemented (IsUnit(lq2)) -- sorry.\n"; return C; elif IsUnit(lq3) then // ??? printf "Not yet implemented (IsUnit(lq3)) -- sorry.\n"; return C; else return Minimise(Corig, BreakUpWith(N, [Integers() | lq2, lq3]) : Check := false); end if; else // split N return Minimise(Corig, BreakUpWith(N, ChangeUniverse(Eltseq(mat), Integers())) : Check := false); end if; // mat eq 0 elif exists{m : m in minors | IsUnit(m)} then // rank is 2 // ??? printf "Not yet implemented (rank(mat) = 2) -- sorry.\n"; return C; else // split N return Minimise(Corig, BreakUpWith(N, ChangeUniverse(minors, Integers())) : Check := false); end if; // all minors = 0 elif not IsUnit(det) then return Minimise(Corig, BreakUpWith(N, [Integers()!det]) : Check := false); end if; // det eq 0 // here mat is invertible mati := mat^-1; pols := [&+[mati[i,j]*derivs[j] : j in [1..3]] : i in [1..3]]; qpols := [pols[i] - PF.1*PF.i : i in [1..3]]; assert forall{q : q in qpols | Derivative(q, 1) eq 0}; p1 := PF.3*qpols[2] - PF.2*qpols[3]; p2 := PF.2*PF.3*qpols[1] + qpols[2]*qpols[3]; if p1 eq 0 and p2 eq 0 then // singular line given by x = -Q2(y,z)/y = -Q3(y,z)/z assert Evaluate(qpols[2], [0, 0, PF.3]) eq 0; Q2 := qpols[2]; return MinLine(C, N, [1, MonomialCoefficient(Q2, PF.2^2), MonomialCoefficient(Q2, PF.2*PF.3)]); end if; gl := [GCD(N, GCD(ChangeUniverse(Coefficients(Fd), Integers()))) : Fd in [p1, p2]]; if exists{g : g in gl | g ne 1} then return Minimise(Corig, BreakUpWith(N, gl) : Check := false); end if; // p1 and p2 are primitive. Find common root. lc1 := MonomialCoefficient(p1, PF.2^3); lc2 := MonomialCoefficient(p2, PF.2^4); if lc1 eq 0 and lc2 eq 0 then // singularity at (-Q3(0,1) : 0 : 1) return MinSing(C, N, [-Evaluate(qpols[3], [0,0,1]), 0, 1]); end if; if not IsUnit(lc1) or not IsUnit(lc2) then return Minimise(Corig, BreakUpWith(N, [Integers() | lc1,lc2]) : Check := false); end if; p1 := 1/lc1*p1; p2 := 1/lc2*p2; ptemp := p2 - PF.2*p1; ptemp := ptemp - MonomialCoefficient(ptemp, PF.2^3*PF.3)*PF.2*p1; p3 := &+[mon*MonomialCoefficient(ptemp, PF.2^2*mon) : mon in [PF.2^i*PF.3^(2-i) : i in [0..2]]]; lc3 := MonomialCoefficient(p3, PF.2^2); if lc3 eq 0 then // lc3 is linear, so must be the GCD lc3 := MonomialCoefficient(p3, PF.2*PF.3); assert lc3 ne 0; if not IsUnit(lc3) then return Minimise(Corig, BreakUpWith(N, [Integers() | lc3]) : Check := false); end if; zsing := -MonomialCoefficient(p3, PF.3^2)/lc3; assert Evaluate(p1, [0,1,zsing]) eq 0 and Evaluate(p2, [0,1,zsing]) eq 0; return MinSing(C, N, [-Evaluate(qpols[2], [0,1,zsing]), 1, zsing]); elif not IsUnit(lc3) then return Minimise(Corig, BreakUpWith(N, [Integers() | lc3]) : Check := false); else p3 := 1/lc3*p3; ptemp := p1 - PF.2*p3; ptemp := ptemp - MonomialCoefficient(ptemp, PF.2^2*PF.3)*PF.2*p3; p4 := &+[mon*MonomialCoefficient(ptemp, PF.2^2*mon) : mon in [PF.2,PF.3]]; lc4 := MonomialCoefficient(p4,PF.2); assert lc4 ne 0; if not IsUnit(lc4) then return Minimise(Corig, BreakUpWith(N, [Integers() | lc4]) : Check := false); end if; zsing := -MonomialCoefficient(p4, PF.3); assert Evaluate(p1, [0,1,zsing]) eq 0 and Evaluate(p2, [0,1,zsing]) eq 0; return MinSing(C, N, [-Evaluate(qpols[2], [0,1,zsing]), 1, zsing]); end if; // lc3 eq 0 end intrinsic; */