//------------------------------------------------------------
// 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 >