misc‎ > ‎

D21

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

Comments