//------------------------------------------------------------ // D21 : Heron triangles with three rational medians //------------------------------------------------------------ P := PolynomialRing(Rationals(),3); Q := PolynomialRing(Rationals(),2); Z := Integers(); h := homQ | [z,1,w]>; // --- parametrization of triangles with two rational medians a := (-2*y*x^2-y^2*x)+(2*x*y-y^2)+x+1; b := (y*x^2+2*y^2*x)+(2*x*y-x^2)-y+1; c := (y*x^2-y^2*x)+(x^2+2*x*y+y^2)+x-y; // --- equation defining the, not necessarily rational, 3rd median M2 := 2*a^2+2*b^2-c^2-4*m^2; // --- 8 curves on which above triangles also have rational area C1 := 27*x^3*y^3-x*y*(x-y)*(8*x^2+11*x*y+8*y^2)-3*x*y*(5*x^2-x*y+5*y^2); C1 := C1-(x-y)*(x^2+4*x*y+y^2)-3*x^2+7*x*y-3*y^2-3*x+3*y-1; C2 := 3*x^2*y^2-2*x*y*(x-y)-x^2-6*x*y-y^2+1; C3 := x*y*(x-y)^3-x^4-11*x^3*y-3*x^2*y^2-11*x*y^3-y^4-2*x^3+2*y^3; C3 := C3+10*x*y+2*x-2*y+1; C4 := x*y*(x-y)+x*y+2*(x-y)-1; C5 := (x-1)^3*y^2+2*(x+1)*(x^3+2*x^2-2*x+1)*y+(2*x-1)*(x+1)^3; C6 := y^4+2*x*y^3-y^3-3*y^2*x-3*y*x^2-2*x*y-x-x^2; C7 := 2*y^4*x-2*y^4+x^2*y^3-6*x*y^3+5*y^3+3*x^2*y^2-3*y^2; C7 := C7+3*y*x^2+2*x*y-y+x^2+2*x+1; C8 := 3*y^2*x-y^2+3*y*x^2-2*x*y+2*x^3*y+y+x^3+x^4; curves := [C1,C2,C3,C4,C5,C6,C7,C8]; for i in [1..8] do // --- genus of intersection V4 : of C=0 and M2=0 eliminating y print "=========================================================";; print "Curve #",i; C := curves[i]; f4 := h(Resultant(M2,C,y)); V4 := Curve(AffineSpace(Q),f4); g := Genus(V4); print V4; print "Genus = ",g; // --- computation of bad primes of V4 f4z := Derivative(f4,z); f4w := Derivative(f4,w); r12 := Resultant(f4,f4z,w); r13 := Resultant(f4,f4w,w); gcd := Gcd(r12,r13); res1 := Resultant(r12 div gcd, r13 div gcd,z); r12 := Resultant(f4,f4z,z); r13 := Resultant(f4,f4w,z); gcd := Gcd(r12,r13); res2 := Resultant(r12 div gcd, r13 div gcd,w); bad_primes := PrimeDivisors(Gcd(Z!res1,Z!res2)); print "Bad primes = ",bad_primes; // ---Weil primes p : p+1-2*g*(Sqrt(p)) <= 0 low := Floor((g-Sqrt(g^2-1))^2); high := Ceiling((g+Sqrt(g^2-1))^2); weil := [p : p in [low..high] | IsPrime(p)]; print "Weil primes = ",weil; // --- find one local solution one_local_solution := function(f4,F) P := PolynomialRing(F,2); f4 := P!f4; for z in F do for w in F do if Evaluate(f4,[z,w]) eq 0 then return {[z,w]}; end if; end for; end for; return {}; end function; // --- find all local solutions all_local_solutions := function(f4,F) P := PolynomialRing(F,2); f4 := P!f4; points := {}; for z in F do for w in F do if Evaluate(f4,[z,w]) eq 0 then if not([z,w] in points) then points := points join {[z,w]}; end if; end if; end for; end for; return points ; end function; // --- look for local solutions for bad primes and Weil primes print "Looking for local solutions ..."; found := 0; S := Sort(Setseq(Seqset(bad_primes cat weil))); for p in S do if #one_local_solution(f4,GF(p)) eq 0 then found := p; end if; end for; if found ne 0 then print "No local solutions at p=",found; else print "Local solutions everywhere."; end if; // --- naive look for affine points lap := function(f4,limit) points := {}; for z_num in [-limit..limit] do for z_den in [1..limit] do for w_num in [-limit..limit] do for w_den in [1..limit] do if Evaluate(f4,[z_num/z_den,w_num/w_den]) eq 0 then points := points join {[z_num/z_den,w_num/w_den]}; end if; end for; end for; end for; end for; return points; end function; // --- faster search for affine points on Aw^4+Bw^2+C=0 sap := function(f4,limit) A := Coefficient(f4,w,4); B := Coefficient(f4,w,2); C := Coefficient(f4,w,0); D := B^2-4*A*C; points := {}; for z_num in [-limit..limit] do for z_den in [1..limit] do if z_num eq 0 or GCD([z_num,z_den]) eq 1 then Z := z_num/z_den; d2 := Rationals()!Evaluate(D,[Z,1]); d2_is_sqr,d := IsSquare(d2); if d2_is_sqr then A2 := Evaluate(2*A,[Z,1]); B2 := Evaluate(B,[Z,1]); if A2 ne 0 then W2 := (d-B2)/A2; W2_is_sqr,W := IsSquare(W2); if W2_is_sqr then points := points join {[Z,W]}; points := points join {[Z,-W]}; end if; W2 := (-d-B2)/A2; W2_is_sqr,W := IsSquare(W2); if W2_is_sqr then points := points join {[Z,W]}; points := points join {[Z,-W]}; end if; else W2 := -Evaluate(C,[Z,1])/B2; W2_is_sqr,W := IsSquare(W2); if W2_is_sqr then points := points join {[Z,W]}; points := points join {[Z,-W]}; end if; end if; end if; end if; end for; end for; return points; end function; if Degree(f4,w) eq 4 and Coefficient(f4,w,3) eq 0 and Coefficient(f4,w,1) eq 0 then print "Looking for rational points fast ..."; points := sap(f4,400); print "Rational points h(z) <= 400"; print points; else print "Looking for rational points slow ..."; points := lap(f4,20); print "Rational points h(z),h(w) <= 20"; print points; end if; // --- Coleman's Bound on V4 : if rank(J_X(Q)) < g, p > 2g and p is good then // #X(Q) = #X(GF(p)) + 2(g-1) // note that [-1,0] is a singular point included in the count // while the point at infinity is excluded so the counts match. good_primes := [p : p in [2*g..100] | not(p in bad_primes) and IsPrime(p)]; min := 10^6; for p in good_primes do solns := #all_local_solutions(f4,GF(p))+2*(g-1); if solns lt min then min := solns; p_min := p; end if; end for; print "Coleman bound = ",min," at p=",p_min; end for; //====================================================================== // --- factor into real quadratics real_factors := function(f) R := Parent(f); roots := Roots(f); norms := Sort([[Norm(roots[i][1]),i] : i in [1..#roots]]); F := []; j := 1; repeat index := Round(norms[j][2]); if IsReal(roots[index][1]) then Append(~F,x-roots[index][1]); j := j+1; else Append(~F,(x-roots[index][1])*(x-roots[index+1][1])); j := j+2; end if; until j gt #roots; return Sort(F); end function; C := ComplexField(); R := PolynomialRing(C); f := x^9+2*x^8-3*x^7+3*x^6-2*x^5+3*x^4-2*x^3+4*x^2-4*x+1; //real_factors(f); |
misc >