Genus2FormalGroup.m010064400007650000024000000345660002675022000143150ustar00gwilliamstaff // Local power series expansions on the Jacobian intrinsic LocalPowerSeries(K::Fld : Precision := 7) -> SeqEnum {} F := FunctionField(K,7 : Global); return LocalPowerSeries([f0,f1,f2,f3,f4,f5,f6] : Precision := Precision); end intrinsic; // RngFrmElt is the dummy type for RngMSerElt intrinsic LocalPowerSeries(ff::[RngElt] : Precision := 7) -> RngFrmElt {} _ , JacRels := JacobianRelations(ff); require #ff eq 7 : "Argument must have length 7."; F := Universe(ff); // We pick up the global object: P := PolynomialRing(F,16 : Global); // We first need to find the 13 relations that express the // s_i in terms of other s_j. These are of the form X00Xij - ... // New polynomial ring "near the origin": R := PolynomialRing(F, 15); // The dehomogenisation homomorphism: dehomog := hom < P -> R | 1, S01, S02, S03, S04, S05, S06, S07, S08, S09, S10, S11, S12, S13, S14, S15>; // We already know S01 and S02. We want to find expressions for the others: sRel := []; sRel[1] := S01; sRel[2] := S02; CurRel := 1; for sIndex in [3..15] do errChk := CurRel; while HomogeneousComponent(dehomog(JacRels[CurRel]), 1) ne R.sIndex do CurRel := CurRel ne #JacRels select (CurRel + 1) else 1; if CurRel eq errChk then // We've cycled through without finding our Si require true: "Error: couldn't find an expression for Si in the relations."; end if; end while; sRel[sIndex] := -HomogeneousComponent(dehomog(JacRels[CurRel]),2); end for; // We now need to recursively generate the power series. // We know that we should reach the required precision in Precision steps. // errChk is to make sure that we're not in an infinite loop. Coef := PolynomialRing(F, 2); R := PowerSeriesRing(Coef); s := [BigO(X^2) : i in [1..15]]; s[1] := X*S1; s[2] := X*S2; errChk := 1; repeat s := [Evaluate(sRel[i], s) : i in [1..15]]; errChk +:= 1; if errChk gt Precision then require true: "Error: couln't reach required precision."; end if; until Minimum([AbsolutePrecision(s[i]) : i in [3..15]]) ge Precision; return s; end intrinsic; intrinsic FindUpperCoordinates(K::Fld : Precision := 8) -> SeqEnum {Returns the last three projective coordinates (ie u13, u14 and u15) when (1:s1:s2:..) is added to (1:t1:t2:..) on the Jacobian corresponding to the curve y^2 = f0 + f1*X + .. + f6*X^6 over K.} F := FunctionField(K,7 : Global); a, b := FindUpperCoordinates([f0,f1,f2,f3,f4,f5,f6] : Prec := Precision); return a, b; end intrinsic; // FIND UPPER COORDINATES // If c = a + b on the Jacobian this intrinsic calculates approximations // for c13, c14 and c15 as power series in K[[s1,s2,t1,t2]]; intrinsic FindUpperCoordinates(ff::[RngElt] : Prec := 8) -> SeqEnum {} // Precision less than 6 gives errors: Prec := Maximum(Prec, 6); // Almost nothing works unless we're over a field!: require IsField(Universe(ff)): "The elements of Argument 1 must belong to a field"; require #ff eq 7: "Argument 1 must have length 7"; // The field that we are working over: F := Universe(ff); f0, f1, f2, f3, f4, f5, f6 := Explode(ff); T := PolynomialRing(F, 30) ; // These functions give the cubic that meets the curve at the two sets of points // given corresponding to (a01,..,a15) and (b01,..,b15). (See Flynn's article) a := []; // mu: mu := s10*t15 + t10*s15 - s14*t11 - t14*s11 + s13*(t12 + t13) + t13*(s12+ s13); // alpha: alpha := s07*t15 + t07*s15 + s13*t09 + t13*s09 - s14*t08 - t14*s08; // beta: beta := -s06*t15 - t06*s15 + s08*(t12 + t13) + t08*(s12 + s13) - s11*t09 - t11*s09; // gamma gamma := s06*t14 + t06*s14 - s07*(t12 + t13) - t07*(s12 + s13) + s09*t10 + t09*s10; // delta delta := -s06*t13 - t06*s13 + s07*t11 + t07*s11 - s08*t10 - t08*s10; // These formulas are lifted straight from Flynn's article: c13 := s15^2*t15^2*(delta^2 - f0*mu^2); c14 := s13*t13*(s15*t15*(f5*mu^2 - 2*alpha*beta) - (s15*t14 + s14*t15)*(alpha^2 - f6*mu^2)); c15 := s13*t13*s15*t15*(alpha^2 - f6*mu^2); PowSer := LocalPowerSeries(ff : Precision:=Prec); Coeffs := PolynomialRing(F,4); P := PowerSeriesRing(Coeffs); Q := Universe(PowSer); B := BaseRing(Q); f := hom< B -> Coeffs | [S1, S2] >; g := hom< B -> Coeffs | [T1, T2] >; s1 := S1*X; s2 := S2*X; t1 := T1*X; t2 := T2*X; aLocalSeries := []; for i in [1..#PowSer] do seq, val := Eltseq(PowSer[i]); p := RelativePrecision(PowSer[i]); aLocalSeries[i] := elt< P | val, [f(seq[j]) : j in [1..#seq]], p>; end for; bLocalSeries := []; for i in [1..#PowSer] do seq, val := Eltseq(PowSer[i]); p := RelativePrecision(PowSer[i]); bLocalSeries[i] := elt< P | val, [g(seq[j]) : j in [1..#seq]], p>; end for; LocalSeries := aLocalSeries cat bLocalSeries; // We now evaluate our exact projective expressions for c13, c14 and c15 // at the series approximations for s1,..,s15 and t1,..,t15: cp13 := Evaluate(c13, LocalSeries); cp14 := Evaluate(c14, LocalSeries); cp15 := Evaluate(c15, LocalSeries); // We take out any common factors to make things nice: seq13 := Eltseq(cp13); seq14 := Eltseq(cp14); seq15 := Eltseq(cp15); DivRenorm := GreatestCommonDivisor(seq13 cat seq14 cat seq15); DivRenorm := DivRenorm*X^(TotalDegree(DivRenorm)); cp13 := P!(cp13/DivRenorm); cp14 := P!(cp14/DivRenorm); cp15 := P!(cp15/DivRenorm); // We know that, after some rescalling we should have: // cp13 ~= (s1+t1)^2*(s2+t2)^2 // cp14 ~= 2*(s1+t1)*(s2+t2)^3 // cp14 ~= (s2+t2)^4 // Here we determine what rescalling is necessary: LeadingTerm := Eltseq(cp13)[1]; ProjectiveFactor := Numerator(LeadingTerm / ((S1+T1)^2*(S2+T2)^2)); MultRenorm := Denominator(LeadingTerm / ((S1+T1)^2*(S2+T2)^2)); MultRenorm := MultRenorm*X^TotalDegree(MultRenorm); ProjectiveFactor := ProjectiveFactor*X^TotalDegree(ProjectiveFactor); cp13 *:= MultRenorm; cp14 *:= MultRenorm; cp15 *:= MultRenorm; return ProjectiveFactor, [cp13, cp14, cp15]; end intrinsic; // PROCEDURE: LinearSolve() // Given a relation, an index, i to solve for (ie if Index = 5 then we // are solving for X05) and the power series expansions so far this // procedure solves for Xi and then clears any denominators. procedure LinearSolve(Rel, Index, ~Exp) // The polynomial ring in X01, X02, X03 etc: M := Parent(Rel); P := Universe(Exp); // We need the power series ring over the function field: Coeffs := BaseRing(P); BiggerCoeffs := FieldOfFractions(Coeffs); Q := PowerSeriesRing(BiggerCoeffs); top, bottom := Explode(Terms(Rel, M.(Index+1))); bottom := ExactQuotient(bottom, M.(Index+1)); print top, bottom; top := Evaluate(top, Exp); bottom := Evaluate(bottom, Exp); print top, bottom; printf "Valuation of top: %o, Valuation of bottom: %o\n", Valuation(top), Valuation(bottom); printf "RelPrec of top: %o ..of bottom: %o\n", RelativePrecision(top), RelativePrecision(bottom); result := Q!(-top / bottom); // We must clear any denominators: seq := Eltseq(result); denom := &*[Denominator(seq[i]) : i in [1..#seq]]; denom := ExactQuotient(denom, GCD([Denominator(seq[i]) : i in [1..#seq]])); denom := denom*X^TotalDegree(denom); // Update all the other coordinates: for i in [1..16] do Exp[i] *:= denom; end for; Exp[Index+1] := P!(denom*result); end procedure; // INTRINSIC: UpdateCoords intrinsic UpdateCoords(ff::[RngElt], AlreadyKnown::[RngIntElt], Wanted::[RngIntElt], ~SeriesExpansions::SeqEnum) {} SolvingSequence := JacobianSolvingSequence(ff, AlreadyKnown, Wanted); // The polynomial ring in X01, X02, X03 etc: M := Parent(SolvingSequence[1][2]); P := Universe(SeriesExpansions); // We need the power series ring over the function field: Coeffs := BaseRing(P); BiggerCoeffs := FieldOfFractions(Coeffs); Q := PowerSeriesRing(BiggerCoeffs); for Step in SolvingSequence do CurrentIndex, Relation, IsSquare := Explode(Step); printf "Attempting to calculate s%o using the following relation:\n%o\n", CurrentIndex, Relation; if not IsSquare then top, bottom := Explode(Terms(Relation, M.(CurrentIndex+1))); bottom := ExactQuotient(bottom, M.(CurrentIndex+1)); top := Evaluate(top, SeriesExpansions); bottom := Evaluate(bottom, SeriesExpansions); result := Q!(-top / bottom); else printf "(Note: We must evaluate a square root)\n"; c, b, a := Explode(Terms(Relation, M.(CurrentIndex+1))); require b eq 0: "Unexpected quadratic! (rather than just a square root)"; a := ExactQuotient(a, (M.(CurrentIndex+1))^2); a := Evaluate(a, SeriesExpansions); c := Evaluate(c, SeriesExpansions); // We need to write a SquareRoot function as this generic // intrinsic doesn't working in characteristic not equal to 0. result := SquareRoot(Q!(-c/a)); end if; // We must clear any denominators: seq := Eltseq(result); denom := &*[Denominator(seq[i]) : i in [1..#seq]]; // We want the lcm: if #seq ne 1 then denom := ExactQuotient(denom, GCD([Denominator(seq[i]) : i in [1..#seq]])); end if; denom := denom*X^TotalDegree(denom); // Update all the other coordinates: for i in {1..16} do SeriesExpansions[i] *:= denom; end for; SeriesExpansions[CurrentIndex+1] := P!(denom*result); printf "Finished (Relative precision: %o. Denominator: %o).\n\n", RelativePrecision(SeriesExpansions[CurrentIndex+1]), Eltseq(denom)[1]; end for; end intrinsic; intrinsic FormalGroup(k::Fld : Precision:=7) -> SeqEnum {} k := FunctionField(k, 7); return FormalGroup([f0,f1,f2,f3,f4,f5,f6] : Precision:= Precision); end intrinsic; intrinsic FormalGroup(ff::[RngElt] : Precision := 7) -> SeqEnum {} _, Jac := JacobianRelations(ff); M := Universe(Jac); f0, f1, f2, f3, f4, f5, f6 := Explode(ff); Expansion := []; // First find x13, x14 and x15: ProjectiveFactor, b := FindUpperCoordinates(ff : Prec := Precision); P := Parent(b[1]); Coeffs := BaseRing(P); s1 := S1*X; s2 := S2*X; t1 := T1*X; t2 := T2*X; BiggerCoeffs := FieldOfFractions(Coeffs); Q := PowerSeriesRing(BiggerCoeffs); Expansion := [ P!1 : i in [1..16]]; Expansion[13 + 1] := b[1]; Expansion[14 + 1] := b[2]; Expansion[15 + 1] := b[3]; // We know that x15/x0 ~= (s2 + t2)^4. // Since x15 ~= (s2 + t2)^4*ProjectiveFactor we know that // x0 ~= ProjectiveFactor. Expansion[0 + 1] := ProjectiveFactor + BigO(X^(Valuation(ProjectiveFactor) + 2)); // We can find x10, x11 and x12 without any need of recursion: UpdateCoords( ff, [13, 14, 15], [12, 11, 10], ~Expansion); printf "Attempting to calculate s5 using a very ugly relation: (This may take some time!)\n"; // The horrible relation that we use to get x5: rel := X05^2*X10*X13 - 1/4*X05^2*X11^2 + f6*X05*X10^2*X13 + 1/2*f5*X05*X10*X11*X13 + f4*X05*X10*X13^2 + 1/2*f3*X05*X11*X13^2 + 1/2*f1*X05*X11*X13*X15 + f2*X05*X13^3 + f0*X05*X13^2*X15 + (f4*f6 - 1/4*f5^2)*X10^2*X13^2 + f3*f6*X10*X11*X13^2 + 1/2*f3*f5*X10*X13^3 - 1/2*f1*f5*X10*X13^2*X15 + f2*f6*X11^2*X13^2 + f1*f6*X11*X12*X13^2 + (f1*f6 + f2*f5)*X11*X13^3 + f0*f6*X12^2*X13^2 + (2*f0*f6 + f1*f5)*X12*X13^3 + f0*f5*X12*X13^2*X14 + (f0*f6 + 2*f1*f5 + f2*f4 - 1/4*f3^2)*X13^4 + (f0*f5 + f1*f4)*X13^3*X14 + 1/2*f1*f3*X13^3*X15 + f0*f4*X13^2*X14^2 + f0*f3*X13^2*X14*X15 + (f0*f2 - 1/4*f1^2)*X13^2*X15^2; c, b, a := Explode(Coefficients(rel, M.6)); printf "Expanding a..."; A := Evaluate(a, Expansion); printf "finished.\nExpanding b..."; B := Evaluate(b, Expansion); printf "finished.\nEvaluation c..."; C := Evaluate(c, Expansion); printf "finished.\nRemoving common factors..."; aseq := Eltseq(A); bseq := Eltseq(B); cseq := Eltseq(C); printf "\n --finished eltseq\n"; gcd := GreatestCommonDivisor(aseq cat bseq cat cseq); printf " --finished gcd\n"; gcd := gcd*X^(TotalDegree(gcd)); A := P!(A/gcd); B := P!(B/gcd); C := P!(C/gcd); printf " --finished quotients\n"; printf "gcd: %o; ...finished\nEvaluating square...", gcd; sq := B^2 - 4*A*C; printf "finished.\nEvaluating root..."; printf "Testing our square root: \n"; return sq, (s2 + t2)^2*Expansion[0+1]; root := OurSquareRoot(Q!sq, (s2 + t2)^2); printf "finished.\n"; printf "Valuation of top: %o, valuation of bottom: %o\n", Valuation(-B + root), Valuation(A); result := Q!((-B - root)/(2*A)); // We must clear any denominators: seq := Eltseq(result); denom := &*[Denominator(seq[i]) : i in [1..#seq]]; // We want the lcm: if #seq ne 1 then denom := ExactQuotient(denom, GCD([Denominator(seq[i]) : i in [1..#seq]])); end if; denom := denom*X^TotalDegree(denom); // Update all the other coordinates: for i in {1..16} do Expansion[i] *:= denom; end for; printf "result*denom = %o\n", result*denom; Expansion[5 + 1] := P!(result*denom); printf "Finished (Relative Precision: %o. Denominator: %o)", RelativePrecision(Expansion[5+1]), denom; UpdateCoords( ff, [15, 14, 13, 12, 11, 10, 5], [0, 9, 8, 7, 6, 4, 3, 2, 1], ~Expansion); return Expansion; end intrinsic; Jacobian.m010064400007650000024000000361250002675022000124730ustar00gwilliamstaff//////////////////////////////////////////////////////////////////////// // Scheme Theoretic Equations // // for Genus 2 Jacobians // //////////////////////////////////////////////////////////////////////// intrinsic JacobianRelations(J::JacHyp) -> SeqEnum {} return JacobianRelations(Curve(J)); end intrinsic; intrinsic JacobianRelations(C::CrvHyp) -> SeqEnum {} f, h := HyperellipticPolynomials(C); require h eq 0 : "Hyperelliptic Curve must be a simplified model. See SimplifiedModel()"; require Genus(C) eq 2 : "Hyperelliptic Curve must have genus 2"; require Degree(f) eq 6 : "At the moment the curve must be defined as y^2 = sextic"; fs := Eltseq(f); return JacobianRelations(fs); end intrinsic; intrinsic JacobianRelations(K::Fld) -> SeqEnum {} F := FunctionField(K,7 : Global); return JacobianRelations([f0,f1,f2,f3,f4,f5,f6]); end intrinsic; intrinsic JacobianRelations(ff::[RngElt]) -> SeqEnum {} require #ff eq 7 : "Argument must have length 7."; F := Universe(ff); f0,f1,f2,f3,f4,f5,f6 := Explode(ff); K := FunctionField(F,4 : Global); t12 := x1 + x2; n12 := x1 * x2; m00 := (y1-y2)/(x1-x2); m01 := (x2*y1-x1*y2)/(x1-x2); m02 := (x2^2*y1-x1^2*y2)/(x1-x2); m03 := (x2^3*y1-x1^3*y2)/(x1-x2); a := [ K | 0 : i in [0..15] ]; a[15+1] := 1; a[14+1] := t12; a[13+1] := n12; a[12+1] := t12^2 - 2*n12; a[11+1] := t12 * n12; a[10+1] := n12^2; a[09+1] := m00; a[08+1] := m01; a[07+1] := m02; a[06+1] := m03; // F0 := 2*f0 + f1*(x1+x2) + 2*f2*(x1*x2) + f3*(x1*x2)*(x1+x2) + 2*f4*(x1*x2)^2 + f5*(x1*x2)^2*(x1+x2) + 2*f6*(x1*x2)^3; // a[05+1] := (F0([x1,x2,y1,y2])-2*y1*y2)/(x1-x2)^2; // F1 := f0*(x1+x2) + 2*f1*(x1*x2) + f2*(x1*x2)*(x1+x2) + 2*f3*(x1*x2)^2 + f4*(x1*x2)^2*(x1+x2) + 2*f5*(x1*x2)^3 + f6*(x1*x2)^3*(x1+x2); a[04+1] := (F1([x1,x2,y1,y2])-(x1+x2)*y1*y2)/(x1-x2)^2; // a[03+1] := a[05+1]*(x1*x2); // G0 := 4*f0 + f1*(x1+3*x2) + 2*f2*(x1+x2)*x2 + f3*(3*x1+x2)*x2^2 + 4*f4*(x1*x2)*x2^2 + f5*(x1*x2)*(x1+3*x2)*x2^2 + 2*f6*(x1*x2)*(x1+x2)*x2^3; a[02+1] := (G0([x1,x2,y1,y2])*y1-G0([x2,x1,y2,y1])*y2)/(x1-x2)^3; // G1 := 2*f0*(x1+x2) + f1*(3*x1+x2)*x2 + 4*f2*(x1*x2)*x2 + f3*(x1*x2)*(x1+3*x2)*x2 + 2*f4*(x1*x2)*(x1+x2)*x2^2 + f5*(x1*x2)*(3*x1+x2)*x2^3 + 4*f6*(x1*x2)^2*x2^3; a[01+1] := (G1([x1,x2,y1,y2])*y1-G1([x2,x1,y2,y1])*y2)/(x1-x2)^3; // a[00+1] := a[05+1]^2; // P := PolynomialRing(F,16 : Global); rels := [ X00*X03 - X01^2 + f4*X03^2 + (-4*f2*f6 + f3*f5)*X03*X10 - 8*f1*f6*X04*X10 - 8*f0*f6*X04*X11 - 4*f0*f5*X04*X13 + f0*X05^2 - f1*f5*X05*X10 + (-4*f1*f5*f6 - 4*f2*f4*f6 + f2*f5^2 + f3^2*f6)*X10^2 + (-4*f0*f5*f6 - 4*f1*f4*f6 + f1*f5^2)*X10*X11 + (-2*f0*f5^2 - 6*f1*f3*f6)*X10*X13 + (-4*f0*f2*f6 - 2*f0*f3*f5 - 3*f1^2*f6)*X10*X15 + (-4*f0*f4*f6 + f0*f5^2)*X11^2 - 8*f0*f3*f6*X11*X13 - 4*f0*f1*f6*X11*X15 - 2*f0*f1*f5*X13*X15, X00*X04 - X01*X02 - f3*f6*X03*X10 + (-9*f0*f5 - f1*f4)*X03*X15 - 4*f2*f6*X04*X10 + (-20*f0*f6 - 3*f1*f5)*X04*X13 + (-7*f1*f6 - f2*f5)*X05*X10 - 2*f0*f4*X05*X14 - f0*f3*X05*X15 - 4*f1*f6*X06*X08 - 4*f0*f6*X06*X09 + 2*f1*f6*X07^2 + 4*f0*f6*X07*X08 - 4*f0*f5*X07*X09 + 4*f0*f5*X08^2 + (-2*f1*f6^2 - 2*f2*f5*f6)*X10^2 + (-4*f0*f6^2 - 2*f1*f5*f6)*X10*X11 - 2*f0*f5*f6*X10*X12 + (-18*f0*f5*f6 - 6*f1*f4*f6 - f1*f5^2 - 2*f2*f3*f6)*X10*X13 + (-8*f0*f4*f6 - f0*f5^2 - 3*f1*f3*f6)*X10*X14 - 4*f0*f3*f6*X12*X13 + (-20*f0*f3*f6 - 4*f0*f4*f5 - 4*f1*f2*f6 - 2*f1*f3*f5)*X13^2 + (-8*f0*f2*f6 - 3*f0*f3*f5 + f1^2*f6)*X13*X14 + (-14*f0*f1*f6 - 6*f0*f2*f5 - f1^2*f5)*X13*X15 + (-4*f0^2*f6 - 2*f0*f1*f5)*X14*X15 - 4*f0^2*f5*X15^2, X00*X05 - X02^2 + f6*X03^2 - f1*f5*X03*X15 - 8*f0*f6*X04*X14 + f2*X05^2 - 4*f0*f5*X05*X14 + (-4*f0*f4 + f1*f3)*X05*X15 - 4*f1*f6*X07*X08 + 4*f1*f6^2*X10*X11 + 2*f1*f5*f6*X10*X13 - 4*f0*f5*f6*X10*X14 + (-4*f0*f4*f6 + f0*f5^2 - 2*f1*f3*f6)*X10*X15 - 8*f0*f3*f6*X11*X15 + (-4*f0*f2*f6 + f1^2*f6)*X12*X15 + (-8*f0*f2*f6 - 2*f0*f3*f5 + 4*f1^2*f6)*X13*X15 + (-4*f0*f2*f5 + f1^2*f5)*X14*X15 + (-4*f0*f2*f4 + f0*f3^2 + f1^2*f4)*X15^2, X00*X06 - X01*X03 - f3*X04*X06 - 2*f2*X05*X06 - 3*f1*X05*X07 - 4*f0*X05*X08 - 4*f2*f6*X06*X10 - 4*f1*f6*X06*X11 - 4*f0*f6*X06*X12 + (-2*f2*f5 + f3*f4)*X07*X10 + (-4*f0*f6 - 4*f1*f5)*X07*X11 - 4*f0*f5*X07*X12 + (4*f0*f6 + 2*f1*f5 + f3^2)*X08*X10 + (-4*f0*f5 - 4*f1*f4 + 3*f2*f3)*X08*X11 - 4*f0*f4*X08*X12 + (-8*f0*f4 - 2*f1*f3 + 4*f2^2)*X08*X13 + (-4*f0*f3 + 2*f1*f2)*X08*X14 + (-4*f0*f2 + 2*f1^2)*X08*X15 + (6*f0*f5 + 4*f1*f4 - 2*f2*f3)*X09*X10 + (4*f0*f4 + f1*f3)*X09*X11 + f0*f3*X09*X12 + (3*f0*f3 + 2*f1*f2)*X09*X13 + 4*f0*f2*X09*X14 + 2*f0*f1*X09*X15, X01*X10 - X03*X06 + f3*X07*X10 + 2*f2*X08*X10 + f1*X08*X11 + f1*X09*X10 + 2*f0*X09*X11, X01*X11 - 2*X04*X06 - f5*X06*X10 + f3*X08*X10 + 2*f2*X08*X11 + f1*X08*X13 + 2*f1*X09*X11 + 2*f0*X09*X12 + 4*f0*X09*X13, X01*X12 - 2*X05*X06 - 4*f6*X06*X10 - f5*X06*X11 - 2*f5*X07*X10 - 2*f4*X07*X11 - f3*X08*X11 + f1*X08*X14 + 2*f0*X09*X14, X01*X13 - X05*X06 + f3*X08*X11 + 2*f2*X08*X13 + f1*X08*X14 - f3*X09*X10 + f1*X09*X13 + 2*f0*X09*X14, X01*X14 - 2*X05*X07 - f5*X07*X11 + f5*X08*X10 - 2*f4*X08*X11 - f3*X08*X13 + f1*X08*X15 + 2*f4*X09*X10 + 2*f0*X09*X15, X01*X15 - X05*X08 + 2*f6*X07*X11 - 2*f6*X08*X10 + f5*X08*X11 - f5*X09*X10, X02*X10 - X04*X06 + f4*X07*X10 + f3*X08*X10 + f2*X08*X11 + f1*X08*X13 + f1*X09*X11 + f0*X09*X12 + 3*f0*X09*X13, X02*X11 - 2*X05*X06 - 2*f6*X06*X10 - f5*X07*X10 + f3*X08*X11 + 2*f2*X08*X13 + 2*f1*X08*X14 - f3*X09*X10 + f1*X09*X13 + 4*f0*X09*X14, X02*X12 - 2*X05*X07 - 2*f6*X06*X11 - 3*f5*X07*X11 + 2*f5*X08*X10 - 4*f4*X08*X11 - f3*X08*X12 - 4*f3*X08*X13 - 2*f2*X08*X14 + 4*f4*X09*X10 + f3*X09*X11 - f1*X09*X14, X02*X13 - X05*X07 + f1*X08*X15 + 2*f0*X09*X15, X02*X14 - 2*X05*X08 + 2*f6*X07*X11 - 2*f6*X08*X10 + f5*X08*X11 - f3*X08*X14 - 2*f2*X08*X15 - f5*X09*X10 + f3*X09*X13 - f1*X09*X15, X03*X07 - X04*X06 + f4*X07*X10 + f3*X08*X10 + f2*X08*X11 + f1*X09*X11 + f0*X09*X12 + f0*X09*X13, X03*X08 - X05*X06 - 2*f6*X06*X10 - f5*X07*X10 + f3*X08*X11 + 2*f2*X08*X13 + f1*X08*X14 - f3*X09*X10 + f1*X09*X13 + 2*f0*X09*X14, X03*X09 - X05*X07 - 2*f6*X06*X11 - 2*f5*X07*X11 + f5*X08*X10 - 2*f4*X08*X11 - f3*X08*X13 + f1*X08*X15 + 2*f4*X09*X10 + 2*f0*X09*X15, X04*X07 - X05*X06 - f6*X06*X10 + f3*X08*X11 + f2*X08*X13 + f1*X08*X14 - f3*X09*X10 + f0*X09*X14, X04*X08 - X05*X07 - f6*X06*X11 - f5*X07*X11 + f5*X08*X10 - f4*X08*X11 + f1*X08*X15 + f4*X09*X10 + f0*X09*X15, X04*X09 - X05*X08 - f6*X06*X12 - f6*X07*X11 - f5*X07*X12 + f6*X08*X10 - f5*X08*X11 - f4*X08*X12 - 2*f4*X08*X13 - f3*X08*X14 - f2*X08*X15 + 2*f5*X09*X10 + f4*X09*X11 + f3*X09*X13, X06*X13 - X07*X11 + X08*X10, X06*X14 - X07*X12 - X08*X11 + 2*X09*X10, X06*X15 - X08*X12 - X08*X13 + X09*X11, X07*X13 - X08*X11 + X09*X10, X07*X14 - X08*X12 - 2*X08*X13 + X09*X11, X07*X15 - X08*X14 + X09*X13, X00*X07 - X02*X03 - f1*X05*X08 - 2*f0*X05*X09, X01*X04 - X02*X03 - f4*X04*X06 - f0*X05*X09 - f3*f6*X06*X10 + (-f3*f5 + f4^2)*X07*X10 - f1*f6*X07*X11 + (f1*f6 - f2*f5 + f3*f4)*X08*X10 + (-f1*f5 + f2*f4)*X08*X11 + (-f0*f5 + f1*f4)*X09*X11 + f0*f4*X09*X12 + f0*f4*X09*X13, X02*X15 - X05*X09 + 2*f6*X07*X12 + 2*f6*X08*X11 + 2*f5*X08*X12 + 3*f5*X08*X13 + 2*f4*X08*X14 + f3*X08*X15 - 4*f6*X09*X10 - 2*f5*X09*X11 - 2*f4*X09*X13, X00*X08 - X02*X04 - f6*X03*X06 - f5*X04*X06 + f2*X05*X08 + f4*f5*X07*X10 + f1*f6*X07*X12 + f3*f5*X08*X10 + (f1*f6 + f2*f5)*X08*X11 + f1*f5*X08*X12 + 2*f1*f5*X08*X13 + f1*f4*X08*X14 + f1*f3*X08*X15 - 2*f1*f6*X09*X10 + f0*f5*X09*X12 + (2*f0*f5 - f1*f4)*X09*X13 + f0*f3*X09*X15, X01*X05 - X02*X04 + f6*X03*X06 + f2*X05*X08 + f1*f6*X07*X12 + f1*f6*X08*X11 + f1*f5*X08*X12 + 2*f1*f5*X08*X13 + f1*f4*X08*X14 + f1*f3*X08*X15 - 2*f1*f6*X09*X10 - f1*f5*X09*X11 + (f0*f5 - f1*f4)*X09*X13 + f0*f3*X09*X15, X00*X09 - X02*X05 - 4*f6*X04*X06 - 3*f5*X05*X06 - 2*f4*X05*X07 - f3*X05*X08 - 4*f5*f6*X06*X10 - f5^2*X07*X10 - 2*f3*f6*X07*X11 + 2*f3*f6*X08*X10 + f3*f5*X08*X11 - 2*f1*f6*X08*X12 + (-2*f1*f6 + 2*f2*f5)*X08*X13 + f1*f5*X08*X14 - f3*f5*X09*X10 + 2*f1*f6*X09*X11 + f1*f5*X09*X13 + 2*f0*f5*X09*X14, X00*X10 - X03^2, X01*X06 - X03^2 - f4*X06^2 - f3*X06*X07 - 2*f2*X07^2 - 3*f1*X07*X08 - 2*f0*X07*X09 - 3*f0*X08^2 + (2*f2*f6 - f3*f5 + f4^2)*X10^2 + (3*f1*f6 - f2*f5 + f3*f4)*X10*X11 + (-2*f0*f6 + 2*f1*f5 - 2*f2*f4 + f3^2)*X10*X13 + (5*f0*f6 - f1*f5 + f2*f4)*X11^2 + (-f0*f5 + f1*f4)*X11*X12 + (3*f0*f5 - f1*f4 + f2*f3)*X11*X13 + f0*f4*X12^2 + f1*f3*X12*X13 + f0*f3*X12*X14 + (-f1*f3 + 2*f2^2)*X13^2 + (-2*f0*f3 + 2*f1*f2)*X13*X14 + (-4*f0*f2 + 2*f1^2)*X13*X15 + 2*f0*f2*X14^2 + 2*f0*f1*X14*X15 + 3*f0^2*X15^2, X03*X10 - X06^2 + f4*X10^2 + f3*X10*X11 + f2*X11^2 + f1*X11*X12 + f1*X11*X13 + f0*X12^2 + 2*f0*X12*X13 + f0*X13^2, X03*X11 - 2*X06*X07 - f5*X10^2 + f3*X10*X13 + 2*f2*X11*X13 + 2*f1*X12*X13 + 2*f0*X12*X14 + 3*f1*X13^2 + 2*f0*X13*X14, X03*X12 - 2*X06*X08 - f5*X10*X11 - 2*f4*X10*X13 - f3*X11*X13 + f1*X13*X14 - 2*f0*X13*X15 + 2*f0*X14^2, X03*X13 - X07^2 + f6*X10^2 + f2*X13^2 + f1*X13*X14 + f0*X14^2, X03*X14 - 2*X07*X08 + 2*f6*X10*X11 + f5*X10*X13 - f3*X13^2 + f1*X13*X15 + 2*f0*X14*X15, X03*X15 - X08^2 + f6*X11^2 + f5*X11*X13 + f4*X13^2 + f0*X15^2, X04*X10 - X06*X07 + f3*X10*X13 + f2*X11*X13 + f1*X12*X13 + f0*X12*X14 + 2*f1*X13^2 + f0*X13*X14, X04*X11 - X06*X08 - X07^2 + f6*X10^2 - f4*X10*X13 + f2*X13^2 + 2*f1*X13*X14 - f0*X13*X15 + 2*f0*X14^2, X04*X12 - X06*X09 - X07*X08 + f6*X10*X11 - f4*X11*X13 - 2*f3*X13^2 - f2*X13*X14 + f0*X14*X15, X04*X13 - X07*X08 + f6*X10*X11 + f5*X10*X13 + f1*X13*X15 + f0*X14*X15, X04*X14 - X07*X09 - X08^2 - f6*X10*X13 + 2*f6*X11^2 + 2*f5*X11*X13 + f4*X13^2 - f2*X13*X15 + f0*X15^2, X05*X10 - X07^2 + f6*X10^2 + f2*X13^2 + f1*X13*X14 + f0*X14^2, X05*X11 - 2*X07*X08 + 2*f6*X10*X11 + f5*X10*X13 - f3*X13^2 + f1*X13*X15 + 2*f0*X14*X15, X05*X12 - 2*X07*X09 - 2*f6*X10*X13 + 2*f6*X11^2 + f5*X11*X13 - f3*X13*X14 - 2*f2*X13*X15 - f1*X14*X15, X05*X13 - X08^2 + f6*X11^2 + f5*X11*X13 + f4*X13^2 + f0*X15^2, X10*X12 + 2*X10*X13 - X11^2, X10*X14 - X11*X13, X10*X15 - X13^2, X11*X14 - X12*X13 - 2*X13^2, X11*X15 - X13*X14, X12*X15 + 2*X13*X15 - X14^2, X00*X11 - 2*X03*X04 - f5*X06^2 - f3*X07^2 - f1*X08^2 + (f3*f6 + f4*f5)*X10^2 + f3*f5*X10*X11 + (f1*f6 + f2*f5)*X11^2 + f1*f5*X11*X12 + 2*f1*f5*X11*X13 + f0*f5*X12^2 + 2*f0*f5*X12*X13 + (f0*f5 + f1*f4 + f2*f3)*X13^2 + f1*f3*X13*X14 + f0*f3*X14^2 + f0*f1*X15^2, X01*X07 - X03*X04 - f1*X08^2 - 2*f0*X08*X09 + f3*f6*X10^2 + 2*f2*f6*X10*X11 + (-f1*f6 + f2*f5)*X10*X13 + 3*f1*f6*X11^2 + 4*f0*f6*X11*X12 + (4*f0*f6 + 2*f1*f5)*X11*X13 + 3*f0*f5*X12*X13 + (5*f0*f5 + f1*f4)*X13^2 + 2*f0*f4*X13*X14 + f0*f3*X13*X15, X02*X06 - X03*X04 - f3*X07^2 - 2*f2*X07*X08 - f1*X07*X09 - 2*f1*X08^2 - 4*f0*X08*X09 + f3*f6*X10^2 + 2*f2*f6*X10*X11 + (-f1*f6 + f2*f5)*X10*X13 + 3*f1*f6*X11^2 + 4*f0*f6*X11*X12 + (4*f0*f6 + 2*f1*f5)*X11*X13 + 3*f0*f5*X12*X13 + (5*f0*f5 + f1*f4)*X13^2 + 2*f0*f4*X13*X14 + f0*f3*X13*X15, X04*X15 - X08*X09 + f6*X11*X12 + f6*X11*X13 + f5*X12*X13 + 2*f5*X13^2 + f4*X13*X14 + f3*X13*X15, X05*X14 - 2*X08*X09 + 2*f6*X11*X12 + 2*f6*X11*X13 + 2*f5*X12*X13 + 3*f5*X13^2 + 2*f4*X13*X14 + f3*X13*X15 - f1*X15^2, X00*X12 - 2*X04^2 - 2*f6*X06^2 - 4*f5*X06*X07 - 2*f4*X07^2 - 4*f3*X07*X08 - 2*f2*X08^2 - 4*f1*X08*X09 - 2*f0*X09^2 + (2*f4*f6 - f5^2)*X10^2 + 4*f3*f6*X10*X11 + 4*f3*f5*X10*X13 + 2*f2*f6*X11^2 + 4*f1*f6*X11*X12 + (4*f1*f6 + 4*f2*f5)*X11*X13 + 2*f0*f6*X12^2 + (4*f0*f6 + 6*f1*f5)*X12*X13 + 4*f0*f5*X12*X14 + (2*f0*f6 + 10*f1*f5 + 2*f2*f4 - f3^2)*X13^2 + (4*f0*f5 + 4*f1*f4)*X13*X14 + 4*f1*f3*X13*X15 + 2*f0*f4*X14^2 + 4*f0*f3*X14*X15 + (2*f0*f2 - f1^2)*X15^2, X00*X13 - X04^2 + f6*X06^2 + f4*X07^2 + f2*X08^2 + f0*X09^2 - f4*f6*X10^2 + f3*f5*X10*X13 - f2*f6*X11^2 - f0*f6*X12^2 + (-2*f0*f6 + f1*f5)*X12*X13 + (-f0*f6 + 2*f1*f5 - f2*f4)*X13^2 + f1*f3*X13*X15 - f0*f4*X14^2 - f0*f2*X15^2, X01*X08 - X04^2 - f6*X06^2 - f5*X06*X07 + f2*X08^2 + f3*f6*X10*X11 + f3*f5*X10*X13 + f2*f6*X11^2 + 2*f1*f6*X11*X12 + (3*f1*f6 + f2*f5)*X11*X13 + 2*f0*f6*X12^2 + (6*f0*f6 + 2*f1*f5)*X12*X13 + 2*f0*f5*X12*X14 + (4*f0*f6 + 4*f1*f5)*X13^2 + (3*f0*f5 + f1*f4)*X13*X14 + f1*f3*X13*X15 + f0*f4*X14^2 + f0*f3*X14*X15, X02*X07 - X04^2 + f4*X07^2 - f1*X08*X09 - f0*X09^2 + f3*f6*X10*X11 + f3*f5*X10*X13 + f2*f6*X11^2 + 2*f1*f6*X11*X12 + (3*f1*f6 + f2*f5)*X11*X13 + 2*f0*f6*X12^2 + (6*f0*f6 + 2*f1*f5)*X12*X13 + 2*f0*f5*X12*X14 + (4*f0*f6 + 4*f1*f5)*X13^2 + (3*f0*f5 + f1*f4)*X13*X14 + f1*f3*X13*X15 + f0*f4*X14^2 + f0*f3*X14*X15, X03*X05 - X04^2 + f6*X06^2 + f4*X07^2 + f2*X08^2 + f0*X09^2 - f4*f6*X10^2 + f3*f5*X10*X13 - f2*f6*X11^2 - f0*f6*X12^2 + (-2*f0*f6 + f1*f5)*X12*X13 + (-f0*f6 + 2*f1*f5 - f2*f4)*X13^2 + f1*f3*X13*X15 - f0*f4*X14^2 - f0*f2*X15^2, X05*X15 - X09^2 + f6*X12^2 + 2*f6*X12*X13 + f5*X12*X14 + f6*X13^2 + f5*X13*X14 + f4*X14^2 + f3*X14*X15 + f2*X15^2, X00*X14 - 2*X04*X05 - f5*X07^2 - f3*X08^2 - f1*X09^2 + f5*f6*X10^2 + f3*f6*X11^2 + f3*f5*X11*X13 + f1*f6*X12^2 + 2*f1*f6*X12*X13 + f1*f5*X12*X14 + (f1*f6 + f2*f5 + f3*f4)*X13^2 + 2*f1*f5*X13*X14 + (f0*f5 + f1*f4)*X14^2 + f1*f3*X14*X15 + (f0*f3 + f1*f2)*X15^2, X01*X09 - X04*X05 - 4*f6*X06*X07 - f5*X06*X08 - 2*f5*X07^2 - 2*f4*X07*X08 - f3*X08^2 + f3*f6*X10*X13 + 2*f2*f6*X11*X13 + 3*f1*f6*X12*X13 + 4*f0*f6*X12*X14 + (5*f1*f6 + f2*f5)*X13^2 + (4*f0*f6 + 2*f1*f5)*X13*X14 + (-f0*f5 + f1*f4)*X13*X15 + 3*f0*f5*X14^2 + 2*f0*f4*X14*X15 + f0*f3*X15^2, X02*X08 - X04*X05 - 2*f6*X06*X07 - f5*X07^2 + f3*f6*X10*X13 + 2*f2*f6*X11*X13 + 3*f1*f6*X12*X13 + 4*f0*f6*X12*X14 + (5*f1*f6 + f2*f5)*X13^2 + (4*f0*f6 + 2*f1*f5)*X13*X14 + (-f0*f5 + f1*f4)*X13*X15 + 3*f0*f5*X14^2 + 2*f0*f4*X14*X15 + f0*f3*X15^2, X00*X15 - X05^2, X02*X09 - X05^2 - 2*f6*X06*X08 - 3*f6*X07^2 - 3*f5*X07*X08 - 2*f4*X08^2 - f3*X08*X09 - f2*X09^2 + 3*f6^2*X10^2 + 2*f5*f6*X10*X11 + (-4*f4*f6 + 2*f5^2)*X10*X13 + 2*f4*f6*X11^2 + f3*f6*X11*X12 + (-2*f3*f6 + 2*f4*f5)*X11*X13 + f2*f6*X12^2 + f3*f5*X12*X13 + (-f1*f6 + f2*f5)*X12*X14 + (-f3*f5 + 2*f4^2)*X13^2 + (3*f1*f6 - f2*f5 + f3*f4)*X13*X14 + (-2*f0*f6 + 2*f1*f5 - 2*f2*f4 + f3^2)*X13*X15 + (5*f0*f6 - f1*f5 + f2*f4)*X14^2 + (3*f0*f5 - f1*f4 + f2*f3)*X14*X15 + (2*f0*f4 - f1*f3 + f2^2)*X15^2 ]; return a, rels; end intrinsic; Relations.m010064400007650000024000000165730002675022000127320ustar00gwilliamstaff//////////////////////////////////////////////////////////////////////// // Verbose mode // //////////////////////////////////////////////////////////////////////// declare verbose RelationsJacG2, 2; //////////////////////////////////////////////////////////////////////// // Generic code for evaluation of polynomials // //////////////////////////////////////////////////////////////////////// intrinsic '@'(S::[RngElt],f::RngMPolElt) -> RngElt {} require #S eq Rank(Parent(f)) : "Argument 1 must have length equal to " * "the rank of the parent of argument 2."; return Evaluate(f,S); end intrinsic; intrinsic '@'(S::[RngElt],f::FldFunRatElt) -> RngElt {} require #S eq Rank(Parent(f)) : "Argument 1 must have length equal to " * "the rank of the parent of argument 2."; return Evaluate(f,S); end intrinsic; // Hack to allow evaluation of rational functions. intrinsic Evaluate(f::FldFunRatMElt, S::[RngElt]) -> RngElt {} return Evaluate(Numerator(f),S) / Evaluate(Denominator(f),S);; end intrinsic; //////////////////////////////////////////////////////////////////////// // Code specific to solving for relations in the genus 2 Jacobian // //////////////////////////////////////////////////////////////////////// intrinsic JacobianWeights() -> SeqEnum {} wts1 := [ 4, 3, 3, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0]; wts2 := [-4,-2,-3, 0,-1,-2, 2, 1, 0,-1, 4, 3, 2, 2, 1, 0]; M := RSpace(Integers(),2); return [ M![wts1[i],wts2[i]] : i in [1..16] ], [ M![2,-i] : i in [0..6] ] cat [ M!0 ]; end intrinsic; intrinsic MonomialsOfWeight(W::[RngIntElt] : Prune := true) -> SeqEnum {} require #W eq 2 : "Argument must be a sequence of two integers."; wts, cfs := JacobianWeights(); w := Universe(wts)!W; wghts := [ [i,j,m,n] : i, j in [1..8], m, n in [1..16] | (cfs[i]+cfs[j]+wts[m]+wts[n]) eq w and i le j and m le n ]; a, r := JacobianRelations(Rationals()); F := Universe(a); K := BaseRing(F); ff := [f0,f1,f2,f3,f4,f5,f6,1]; P := Universe(r); if not Prune then return [ ff[w[1]]*ff[w[2]]*a[w[3]]*a[w[4]] : w in wghts ], [ ff[w[1]]*ff[w[2]]*P.w[3]*P.w[4] : w in wghts ]; end if; ainvs := [ ]; terms := [ ]; mons := {@ @}; for w in wghts do x, y, i, j := Explode(w); m := P.i*P.j; if m notin mons then Append(~ainvs,ff[x]*ff[y]*a[i]*a[j]); Append(~terms,ff[x]*ff[y]*m); Include(~mons,m); end if; end for; return ainvs, terms; end intrinsic; intrinsic MonomialWeight(m::RngMPolElt) -> SeqEnum {} I := [ ]; E := Exponents(m); for i in [1..16] do if E[i] ge 1 then Append(~I,i); end if; if E[i] ge 2 then Append(~I,i); end if; end for; i, j := Explode(I); wts, cfs := JacobianWeights(); wi1 := wts[i][1]; wj1 := wts[j][1]; wi2 := wts[i][2]; wj2 := wts[j][2]; c := Coefficients(m); assert #c eq 1; E := Exponents(Monomials(Numerator(c[1]))[1]); F := Exponents(Monomials(Denominator(c[1]))[1]); c1 := 2*&+E - 2*&+F; c2 := &+[ -(i-1)*(E[i]-F[i]) : i in [1..7] ]; vprintf RelationsJacG2, 2 : "W = [%2o,%2o] + [%2o,%2o] + [%2o,%2o]\n", wts[i][1], wts[i][2], wts[j][1], wts[j][2], c1, c2; W := [ wi1+wj1+c1, wi2+wj2+c2 ]; return W; end intrinsic; intrinsic MonomialRelationsOfWeight(W::[RngIntElt] : Relations := [], StrictWeight := true) -> SetIndx {} require #W eq 2 : "Argument must have length 2."; ainvs, mons := MonomialsOfWeight(W); V := LinearRelations(ainvs : Relations := Relations); WghtRels := {@ @}; wts, cfs := JacobianWeights(); for v in Basis(V) do f := &+[ v[i]*mons[i] : i in [1..#mons] ]; k := Random([ i : i in [1..Degree(V)] | v[i] ne 0 ]); wghtf := MonomialWeight(v[k]*mons[k]); if not StrictWeight or W eq wghtf then Include(~WghtRels,f); else vprintf RelationsJacG2 : "Found relation of weight [%2o,%2o] (not [%2o,%2o]):\n", wghtf[1], wghtf[2], W[1], W[2]; vprint RelationsJacG2 : f; end if; end for; return WghtRels; end intrinsic; /* F := FunctionField(Rationals(),7 : Global); P := PolynomialRing(F,16 : Global); a, r := JacobianRelations([f0,f1,f2,f3,f4,f5,f6]); R := PolynomialRing(F,4 : Global); F1 := y1^2-(f6*x1^6 + f5*x1^5 + f4*x1^4 + f3*x1^3 + f2*x1^2 + f1*x1 + f0); F2 := y2^2-(f6*x2^6 + f5*x2^5 + f4*x2^4 + f3*x2^3 + f2*x2^2 + f1*x2 + f0); DomRels := [F1,F2]; I := ideal< R | DomRels >; //time MonomialRelationsOfWeight([4,0] : Relations := DomRels); // TotRels := {@ @}; for i in [6..0 by -1] do if i eq 6 then jInt := [-4..-6 by -1]; elif i eq 5 then jInt := [-2..-5 by -1]; elif i eq 4 then jInt := [ 0..-4 by -1]; elif i eq 3 then jInt := [ 2..-3 by -1]; elif i eq 2 then jInt := [ 4..-2 by -1]; elif i eq 1 then jInt := [ 3.. 1 by -1]; elif i eq 0 then jInt := [ 6.. 2 by -1]; end if; for j in jInt do printf "Weight: [%2o,%2o]\n", i, j; NewRels := MonomialRelationsOfWeight([i,j] : Relations := DomRels); print NewRels; TotRels join:= {@ Factorization(f)[1][1] : f in NewRels @}; print "Total relations:", #TotRels; end for; end for; */ //////////////////////////////////////////////////////////////////////// // Completely generic... //////////////////////////////////////////////////////////////////////// intrinsic AlgebraicRelations(S::[FldFunRatElt] : DomainRelations := [], CodomainRelations := []) -> RngMPol {} m := #S; R := Universe(S); n := Rank(R); P := PolynomialRing(R,m+n); AssignNames(~P, [ "r" cat IntegerToString(i): i in [1..m]] cat [ "x" cat IntegerToString(i): i in [1..n]]); vars := [ P.j : j in [m+1..m+n] ]; rels := [ Denominator(S[i])(vars)*P.i - Numerator(S[i])(vars) : i in [1..m] ]; if #DomainRelations ne 0 then Q := Universe(DomainRelations); require Type(Q) eq RngMPol : "Parameter DomainRelations must consist of polynomials."; require Rank(Q) eq n : "Rank of universe of DomainRelations " * "must equal universe of argument 1."; rels cat:= [ f([ P.j : j in [m+1..m+n] ]) : f in Relations ]; end if; return EliminationIdeal(ideal< P | rels >,{P.i : i in [1..m]}); end intrinsic; intrinsic LinearRelations(S::[FldFunRatElt] : Relations := []) -> ModTupFld {} for i in [1..#S] do den := Denominator(S[i]); if den ne 1 then S := [ den*f : f in S ]; end if; end for; S := [ Numerator(f) : f in S ]; if #Relations ne 0 then vprint RelationsJacG2 : "Putting functions in normal form."; P := Universe(Relations); require P eq Universe(S) : "Invalid Relations parameter."; I := ideal< Universe(S) | Relations >; S := [ NormalForm(f,I) : f in S ]; end if; vprint RelationsJacG2 : "Extracting coefficients for matrix."; K := BaseRing(Universe(S)); mons := &join[ SequenceToSet(Monomials(f)) : f in S ]; vprint RelationsJacG2 : "Number of monomials:", #mons; vprintf RelationsJacG2 : "Matrix size %o x %o\n", #S, #mons; M := Matrix([ [ MonomialCoefficient(f,m) : m in mons ] : f in S ]); if GetVerbose("RelationsJacG2") ge 2 then for j in [1..#mons] do "j =", j; for i in [1..#S] do if M[i,j] ne 0 then printf "M[%3o,%2o] = %o\n", i, j, M[i,j]; end if; end for; end for; end if; return Kernel(M); end intrinsic; RelationsRoutines.m010064400007650000024000000157130002675022000144560ustar00gwilliamstaff intrinsic FindAllRels(ff::[RngElt], NecessaryMonomials, PossibleMonomials) -> SeqEnum {} _, JacRels := JacobianRelations(ff); M := Universe(JacRels); F := BaseRing(M); // We need to apply vector space routines: require IsField(F): "Base ring must be a field"; MonList := [M|]; for i in [1..16] do MonList cat:= [ M.i*M.j : j in [i..16]]; end for; V := VectorSpace(F, #MonList); RelsMat := []; for i in [1..#JacRels] do Mons := Monomials(JacRels[i]); Coeffs := Coefficients(JacRels[i]); RelsMat[i] := &+[Coeffs[j]*V.(Index(MonList,Mons[j])) : j in [1..#Mons]]; end for; RelsMat := Matrix(RelsMat); // We first list all the monomials that can appear in our relation: WantedMonsVectors := [ V.Index(MonList, mon) : mon in PossibleMonomials]; // We then create the corresponding space... WantedSpace := sub < V | WantedMonsVectors >; // ... and intersect it with the relations: WantedSpace := WantedSpace meet Image(RelsMat); // We now list all the monomials (one of which _must_ appear); NeededMons := [ V.Index(MonList, mon) : mon in NecessaryMonomials]; B := Basis(WantedSpace); PossibleRelations := {}; for i in [1..#B] do for j in NeededMons do if InnerProduct( j , B[i]) ne 0 then sol := Solution(RelsMat, B[i]); sol := Eltseq(sol); PossibleRelations join:= {&+[sol[i]*JacRels[i] : i in [1..#JacRels]]}; end if; end for; end for; // Have we found a relation? if not IsEmpty(PossibleRelations) then return SetToSequence(PossibleRelations); else return false; end if; end intrinsic; // JACOBIAN SOLVING RELATIONS intrinsic JacSolveRels(ff::[RngElt], AlreadyKnown, Wanted) -> SeqEnum {} _, JacRels := JacobianRelations(ff); M := Universe(JacRels); F := BaseRing(M); require IsField(F): "Base ring must be a field"; MonList := [M|]; for i in [1..16] do MonList cat:= [ M.i*M.j : j in [i..16]]; end for; V := VectorSpace(F, #MonList); RelsMat := []; for i in [1..#JacRels] do Mons := Monomials(JacRels[i]); Coeffs := Coefficients(JacRels[i]); RelsMat[i] := &+[Coeffs[j]*V.(Index(MonList,Mons[j])) : j in [1..#Mons]]; end for; RelsMat := Matrix(RelsMat); // We first list all the monomials that can appear in our relation: WantedMons := SetToSequence({ M.(i+1)*M.(j+1) : i in AlreadyKnown, j in AlreadyKnown}); WantedMons cat:= [ M.(i+1)*M.(Wanted+1) : i in AlreadyKnown]; WantedMons cat:= [ M.(i+1)*M.(Wanted + 1) : i in AlreadyKnown]; WantedMonsVectors := [ V.Index(MonList, mon) : mon in WantedMons]; // We then create the corresponding space... WantedSpace := sub < V | WantedMonsVectors >; // ... and intersect it with the relations: WantedSpace := WantedSpace meet Image(RelsMat); // We now need to find an element in the WantedSpace that // contains a term of the form S.CurrentIndex*S.j for some j > CurrentIndex. NeededMons := [ M.(Wanted+1)*M.(i+1) : i in AlreadyKnown]; B := Basis(WantedSpace); PossibleRelations := []; for i in [1..#B] do for j in NeededMons do if InnerProduct( V.Index(MonList, j), B[i]) ne 0 then sol := Solution(RelsMat, B[i]); sol := Eltseq(sol); PossibleRelations cat:= [&+[sol[i]*JacRels[i] : i in [1..#JacRels]]]; end if; end for; end for; if not IsEmpty(PossibleRelations) then // We've found a relation: // "false" because we don't need to evaluate a square root to solve: return PossibleRelations, false; end if; // We are here if, regrettably(!), we need to evaluate a square root: // We add the square to both wanted and needed: WantedMonsVectors cat:= [ V.Index(MonList, M.(Wanted+1)*M.(Wanted+1))]; NeededMons cat:= [ M.(Wanted+1)*M.(Wanted+1)]; WantedSpace := Image(RelsMat) meet sub < V | WantedMonsVectors >; B := Basis(WantedSpace); PossibleRelations := []; for i in [1..#B] do for j in NeededMons do if InnerProduct( V.Index(MonList, j), B[i]) ne 0 then sol := Solution(RelsMat, B[i]); sol := Eltseq(sol); PossibleRelations cat:= [&+[sol[i]*JacRels[i] : i in [1..#JacRels]]]; end if; end for; end for; if not IsEmpty(PossibleRelations) then return PossibleRelations, true; else printf "Failed to find a relation for s%o\n", Wanted; return false; end if; end intrinsic; // JACOBIAN SOLVING SEQUENCE intrinsic JacobianSolvingSequence(ff::[RngElt], AlreadyKnown, Wanted) -> SeqEnum {} ReturnSequence := []; // We return a sequence in which each entry has the form: // < index, relation, true/false>; // Where index is the s_i that it is possible to solve for, relation is the // relation from which it is possible to derive s_i and true/false // indicates whether or not it is necessary to evaluate a square root: AK := AlreadyKnown; for x in [1..#Wanted] do CurrentIndex := Wanted[x]; Rels, NeedSquareRoot := JacSolveRels(ff, AK, CurrentIndex); Sizes := [ 1/#Coefficients(Rels[i]) : i in [1..#Rels]]; _, RelIndex := Maximum(Sizes); ReturnSequence[x] := ; AK cat:= [Wanted[x]]; end for; return ReturnSequence; end intrinsic; intrinsic x5Rel(ff::[RngElt]) -> SeqEnum {} _, Jac := JacobianRelations(ff); M := Universe(Jac); N := FieldOfFractions(M); Expand := [N.i : i in [1..16]]; rels := JacSolveRels(ff, [15,14,13,12,11,10,5], 0); pos0 := []; for i in [1..#rels] do top, bottom := Explode(Coefficients(rels[i],M.1)); pos0[i] := -(top/bottom); end for; rels := JacSolveRels(ff, [15,14,13,12,11,10,5], 3); pos3 := []; for i in [1..#rels] do top, bottom := Explode(Coefficients(rels[i],M.4)); pos3[i] := -(top/bottom); end for; rels := JacSolveRels(ff, [15,14,13,12,11,10,5], 4); pos4 := []; for i in [1..#rels] do top, bottom := Explode(Coefficients(rels[i], M.5)); pos4[i] := -(top/bottom); end for; Rels := {}; rels5 := JacSolveRels(ff, [15, 14, 13, 12, 11, 10, 3, 4], 5); for rel0 in pos0, rel3 in pos3, rel4 in pos4, rel in rels5 do Expand[1] := rel0; Expand[4] := rel3; Expand[5] := rel4; Rels join:= {Numerator(Evaluate(rel, Expand))}; end for; return Rels; end intrinsic;