//---------------------------------------------------
// D21 : Heron triangles with two rational medians
// Ralph Buchholz
//---------------------------------------------------

Z := Integers();
doh := procedure()
  print "ABKLT(P) --> generates orbit of P";
end procedure;

// new Dual Hayes parameters
new_delta := function(theta,phi)
  Theta := (theta*phi+2*phi^2-theta-phi-1)/(3*theta*phi+theta-phi+1);
  Phi   := (-2*theta^2-theta*phi-theta-phi+1)/(3*theta*phi+theta-phi+1);
  return Theta,Phi;
end function;

// negative Hayes parameters
nu := function(theta,phi)
  return -phi,-theta;
end function;

// half dual Hayes parameters
chi := function(theta,phi)
  Theta := (theta*phi+2*phi^2-theta-phi-1)/(3*theta*phi+theta-phi+1);
  return Theta,phi;
end function;

// half dual Hayes parameters
chi_theta := function(theta,phi)
  Theta := (2*theta^2+theta*phi+theta+phi-1)/(3*theta*phi+theta-phi+1);
  return Theta,-theta;
end function;

// half dual Hayes parameters
chi_phi := function(theta,phi)
  Phi := (-theta*phi-2*phi^2+theta+phi+1)/(3*theta*phi+theta-phi+1);
  return -phi,Phi;
end function;

// Old Dual Hayes parameters
delta := function(theta,phi)
  Theta := (2*theta^2+theta*phi+theta+phi-1)/(3*theta*phi+theta-phi+1);
  Phi   := (-theta*phi-2*phi^2+theta+phi+1)/(3*theta*phi+theta-phi+1);
  return Theta,Phi;
end function;

// interchange the sign of a
signa := function(theta,phi)
  Theta := (1-theta)/(1+theta);
  Phi   := -1/phi;
  return Theta,Phi;
end function;

// interchange the sign of b
signb := function(theta,phi)
  Theta := -1/theta;
  Phi   := (phi+1)/(phi-1);
  return Theta,Phi;
end function;

// generate all Hayes parameters
hayes_group := function(theta,phi)
  t1,p1 := Explode([theta,phi]);
  t2,p2 := delta(t1,p1);
  t3,p3 := chi(t1,p1);
  t4,p4 := delta(t3,p3);
  t5,p5 := nu(t1,p1);
  t6,p6 := nu(t2,p2);
  t7,p7 := nu(t3,p3);
  t8,p8 := nu(t4,p4);
  orbit1:= [[t1,p1],[t2,p2],[t3,p3],[t4,p4],[t5,p5],[t6,p6],[t7,p7],[t8,p8]];

  t9,p9 := signa(t1,p1);
  tA,pA := delta(t9,p9);
  tB,pB := chi(t9,p9);
  tC,pC := delta(tB,pB);
  tD,pD := nu(t9,p9);
  tE,pE := nu(tA,pA);
  tF,pF := nu(tB,pB);
  tG,pG := nu(tC,pC);
  orbit2:= [[t9,p9],[tA,pA],[tB,pB],[tC,pC],[tD,pD],[tE,pE],[tF,pF],[tG,pG]];
  
  tH,pH := signb(t1,p1);
  tI,pI := delta(tH,pH);
  tJ,pJ := chi(tH,pH);
  tK,pK := delta(tJ,pJ);
  tL,pL := nu(tH,pH);
  tM,pM := nu(tI,pI);
  tN,pN := nu(tJ,pJ);
  tO,pO := nu(tK,pK);
  orbit3:= [[tH,pH],[tI,pI],[tJ,pJ],[tK,pK],[tL,pL],[tM,pM],[tN,pN],[tO,pO]];

  tP,pP := signb(t9,p9);
  tQ,pQ := delta(tP,pP);
  tR,pR := chi(tP,pP);
  tS,pS := delta(tR,pR);
  tT,pT := nu(tP,pP);
  tU,pU := nu(tQ,pQ);
  tV,pV := nu(tR,pR);
  tW,pW := nu(tS,pS);
  orbit4:= [[tP,pP],[tQ,pQ],[tR,pR],[tS,pS],[tT,pT],[tU,pU],[tV,pV],[tW,pW]];

  return orbit1 cat orbit2 cat orbit3 cat orbit4;
end function;

// Hayes parameters to sides
hayes2sides := function(theta,phi)
  a := (-2*phi*theta^2-phi^2*theta)+(2*theta*phi-phi^2)+theta+1;
  b := (phi*theta^2+2*phi^2*theta)+(2*theta*phi-theta^2)-phi+1;
  c := (phi*theta^2-phi^2*theta)+(theta^2+2*theta*phi+phi^2)+theta-phi;
  A := Denominator(a);
  B := Denominator(b);
  C := Denominator(c);
  l := Lcm([A,B,C]);
  a,b,c := Explode([Z!(a*l),Z!(b*l),Z!(c*l)]);
  g := Gcd([a,b,c]);
  return [a div g,b div g,c div g];
end function;


// Sides to Hayes parameters
sides2hayes := function(a,b,c)
  k := 2*a^2+2*c^2-b^2;
  l := 2*b^2+2*c^2-a^2;
  if (IsSquare(k)) then
    t1 := (c-a+Isqrt(k))/(a+b+c);
    t2 := (c-a-Isqrt(k))/(a+b+c);
  else
    t1 := (c-a+Sqrt(k))/(a+b+c);
    t2 := (c-a-Sqrt(k))/(a+b+c);
  end if;
  if (IsSquare(l)) then
    p1 := (b-c+Isqrt(l))/(a+b+c);
    p2 := (b-c-Isqrt(l))/(a+b+c);
  else
    p1 := (b-c+Sqrt(l))/(a+b+c);
    p2 := (b-c-Sqrt(l))/(a+b+c);
  end if;
  return [t1,p1],[t2,p2],[t1,p2],[t2,p1];
end function;


area_septic := function(theta,phi)
  P := PolynomialRing(Rationals(),2);
  gamma := x*y*(1-x^2)*(1-y^2)*(3*x*y+x-y+1)*(2*x+y-1)*(x+2*y+1)*(x-y+1);
  g := Evaluate(gamma,[theta,phi]);
  return g;
end function;
  

A := function(P)
  theta,phi := Explode(P);
  Ax := -(-theta*phi-theta+phi^2-1)/(2*theta-1+phi)/phi;
  Ay := (3*theta*phi+theta+1-phi)/(2*theta^2+theta*phi+theta+phi-1);
  return ;
end function;
B := function(P)
  theta,phi := Explode(P);
  Bx := -(3*theta*phi+theta+1-phi)/(theta*phi+2*phi^2-phi-1-theta);
  By := -(-phi+theta*phi-theta^2+1)/(2*phi+theta+1)/theta;
  return ;
end function;
C := function(P)
  theta,phi := Explode(P);
  Cx := (2*theta-1+phi)*phi/(-theta*phi-theta+phi^2-1);
  Cy := (2*phi+theta+1)*theta/(-phi+theta*phi-theta^2+1);
  return ;
end function;
K := function(P)
  theta,phi := Explode(P);
  Kx := theta;
  Ky := -(2*theta^2+theta*phi+theta+phi-1)/(3*theta*phi+theta+1-phi);
  return ;
end function;
L := function(P)
  theta,phi := Explode(P);
  Lx := (theta*phi+2*phi^2-phi-1-theta)/(3*theta*phi+theta+1-phi);
  Ly := phi;
  return ;
end function;
T := function(P)
  theta,phi := Explode(P);
  Tx := (2*theta^2+theta*phi+theta+phi-1)/(3*theta*phi+theta+1-phi);
  Ty := -(theta*phi+2*phi^2-phi-1-theta)/(3*theta*phi+theta+1-phi);
  return ;
end function;

ABKLT := function(P)
  image := [];
  Append(~image,P);
  Append(~image,T(P));
  Append(~image,L(P));
  Append(~image,T(L(P)));

  Append(~image,K(P));
  Append(~image,T(K(P)));
  Append(~image,L(K(P)));
  Append(~image,T(L(K(P))));

  Append(~image,B(P));
  Append(~image,T(B(P)));
  Append(~image,L(B(P)));
  Append(~image,T(L(B(P))));

  Append(~image,K(B(P)));
  Append(~image,T(K(B(P))));
  Append(~image,L(K(B(P))));
  Append(~image,T(L(K(B(P)))));

  Append(~image,A(P));
  Append(~image,T(A(P)));
  Append(~image,L(A(P)));
  Append(~image,T(L(A(P))));

  Append(~image,K(A(P)));
  Append(~image,T(K(A(P))));
  Append(~image,L(K(A(P))));
  Append(~image,T(L(K(A(P)))));

  Append(~image,B(A(P)));
  Append(~image,T(B(A(P))));
  Append(~image,L(B(A(P))));
  Append(~image,T(L(B(A(P)))));

  Append(~image,K(B(A(P))));
  Append(~image,T(K(B(A(P)))));
  Append(~image,L(K(B(A(P)))));
  Append(~image,T(L(K(B(A(P))))));
  return image;
end function;
  
lines := function(points)
  Q := [1,-1];
  lines := [];
  for P in points do
    m := (P[2]-Q[2])/(P[1]-Q[1]);
    b := Q[2]-m*Q[1];
    Append(~lines,);
  end for;
  return lines;
end function;

smallest := function(list)
  min := Max(Abs(Numerator(list[1,2])),Abs(Denominator(list[1,2])));
  P_min := list[1];
  for P in list do
    t := Max(Abs(Numerator(P[2])),Abs(Denominator(P[2])));
    if (t lt min and Abs(P[2]) lt 1) then
      P_min := P;
      min := t;
    end if;
  end for;
  return P_min;
end function;

C1 := function(x,y)
  f := 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);
  f := f-(x-y)*(x^2+4*x*y+y^2)-3*x^2+7*x*y-3*y^2-3*x+3*y-1;
  if (f eq 0) then
    return true;
  else
    return false;
  end if;
end function;

C2 := function(x,y)
  f := 3*x^2*y^2-2*x*y*(x-y)-x^2-6*x*y-y^2+1;
  if (f eq 0) then
    return true;
  else
    return false;
  end if;
end function;

C3 := function(x,y)
  f := 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;
  f := f+10*x*y+2*x-2*y+1;
  if (f eq 0) then
    return true;
  else
    return false;
  end if;
end function;

C4 := function(t,f)
  if (t*f*(t-f)+t*f+2*(t-f)-1 eq 0) then
    return true;
  else
    return false;
  end if;
end function;

C5 := function(x,y)
  f := (x-1)^3*y^2+2*(x+1)*(x^3+2*x^2-2*x+1)*y+(2*x-1)*(x+1)^3;
  if (f eq 0) then
    return true;
  else
    return false;
  end if;
end function;

C6 := function(x,y)
  f := y^4+2*x*y^3-y^3-3*y^2*x-3*y*x^2-2*x*y-x-x^2;
  if (f eq 0) then
    return true;
  else
    return false;
  end if;
end function;

C7 := function(x,y)
  f := 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;
  f := f+3*y*x^2+2*x*y-y+x^2+2*x+1;
  if (f eq 0) then
    return true;
  else
    return false;
  end if;
end function;

C8 := function(x,y)
  f := 3*y^2*x-y^2+3*y*x^2-2*x*y+2*x^3*y+y+x^3+x^4;
  if (f eq 0) then
    return true;
  else
    return false;
  end if;
end function;

    
is_sporadic := function(theta,phi)
  orbit := ABKLT();
  for P in orbit do
    x,y := Explode(P);
    if (C1(x,y) or C2(x,y) or C3(x,y) or C4(x,y) or
        C5(x,y) or C6(x,y) or C7(x,y) or C8(x,y)) then
      return false;
    end if;
  end for;
  return true;
end function;


//--------------------------------------------------------------------
// Eulerian 3 median parametrization 
eulerian_3med := function(x)
  a := (2*x*(-9*x^4+10*x^2+3));
  b := ((9*x^4-6*x^2+1)-x*(9*x^4+26*x^2+1));
  c := ((9*x^4-6*x^2+1)+x*(9*x^4+26*x^2+1));
  p := (a+b+c);
  A2 := p*(p-2*a)*(p-2*b)*(p-2*c);
  return A2;
end function;


// --- create naive curve
print "\n------Model------";
P   := PolynomialRing(Integers());
Q := PolynomialRing(Rationals(),2);
f := eulerian_3med(x) div x^2;
h := hom
 Q | u>;
C := Curve(AffineSpace(Q),v^2-h(f));
print "Genus=",Genus(C)," degree=",Degree(f);


// --- use covering map u^2 --> u
c := Coefficients(f);
g := &+[c[2*i-1]*x^(i-1) : i in [1..(#c+1) div 2] ];
C := Curve(AffineSpace(Q),v^2-h(g));
print "Genus=",Genus(C)," degree=",Degree(g);


// --- lower degree by 1
FF := FunctionField(Rationals());
g2 := g div 64;
G2 := FF!g2;
g3 := P!Numerator(Evaluate(G2,1+1/(X)));
C  := Curve(AffineSpace(Q),v^2-h(g3));
print "Genus=",Genus(C)," degree=",Degree(g3);
H3 := HyperellipticCurve(g3);
J3 := Jacobian(H3);
print "J(Q)=",J3;


// --- find rational points (via Stoll only found 7 up to height 10^5)

P0 := [1,0,0];
P1 := [ -9 / 2, 2^5*3^6 ];
P2 := [ -3 / 2, 2^3*3^2 ];
P3 := [ -1 / 2, 0 ];
P4 := [ -7 / 6, 2^5/3^2 ];
P5 := [-9 / 8, 0 ];
P6 := [-9 / 10, 0 ];
points := [H3!P1,H3!P2,H3!P3,H3!P4,H3!P5,H3!P6];
print "points=",points;


// --- build divisors in jacobian from points
D0 := J3![H3!P0,H3!P1];  // gen
D1 := J3![H3!P0,H3!P2];  // gen
D2 := J3![H3!P0,H3!P3];  // gen order 2
D3 := J3![H3!P0,H3!P4];  // D0+4*D1+D3 --> D21 order 2
D4 := J3![H3!P0,H3!P5];  // gen order 2
D5 := J3![H3!P0,H3!P6];  // gen order 2
D6 := J3![H3!P1,H3!P2];  // -D0+D1
D7 := J3![H3!P1,H3!P3];  // -D0-D2
D8 := J3![H3!P1,H3!P4];  // -D0+D3
D9 := J3![H3!P1,H3!P5];  // -D0+D4;
D10 := J3![H3!P1,H3!P6]; // -D0+D5;
D11 := J3![H3!P2,H3!P3]; // -D1+D2
D12 := J3![H3!P2,H3!P4]; // -D1+D3
D13 := J3![H3!P2,H3!P5]; // -D1+D4
D14 := J3![H3!P2,H3!P6]; // -D1+D5
D15 := J3![H3!P3,H3!P4]; // D3+D2
D16 := J3![H3!P3,H3!P5]; // D2+D4; order 2
D17 := J3![H3!P3,H3!P6]; // D2+D5; order 2
D18 := J3![H3!P4,H3!P5]; // -D3+D4
D19 := J3![H3!P4,H3!P6]; // -D3+D5
D20 := J3![H3!P5,H3!P6]; // D4+D5; order 2
D21 := D0+4*D1+D3;       // gen order 2

// build divisor from factorization of g3
D22 :=J3![12*x^2 + 4*x - 9,0]; // D2+D4+D5+D21; order 2

gens := [D0,D1,D2,D4,D5,D21,D22];
D := [D0,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13,D14,D15,D16,D17,D18,D19,D20,D21,D22];
print "Divisors=",gens;


// --- find the bad primes
bad_primes := Factorization(Integers()!Discriminant(H3));


// --- look for torsion subgroup of jacobian (seems to be <= 64)
print "\n------Torsion------";
T := 0;
for p in [i : i in [11..41]|IsPrime(i)] do
  P11 := PolynomialRing(GF(p));
  h_U := hom
 P11 | U>;
  H   := HyperellipticCurve(h_U(g3));
  J   := Jacobian(H);
  T   := GCD(T,#J);
  print p,"Torsion <=",T;
end for;
print "Torsion <=",T;


// --- map the 16 order 2 points into C_2 subgroups of J(GF(p)) for p > 7
// --- a la ZFB ... (this forces #J_tors(Q) <= 16)
SetSeed(76543);
found := [ : i in [1..16]];
for p in [i : i in [11..31]|IsPrime(i)] do
  P_p := PolynomialRing(GF(p));
  h_p := hom
 P_p | U>;
  H_p := HyperellipticCurve(h_p(g3));
  J_p := Jacobian(H_p);
  G,m := AbelianGroup(J_p);
  G_2 := [ m(g) : g in Generators(G) | Order(g) eq 2 ];
  H_2 := [];
  for b in [0..2^#G_2-1] do
    bits := IntegerToSequence(b,2);
    bits := bits cat [0 : i in [1..#G_2-#bits]];
    elt  := &+[bits[i]*G_2[i] : i in [1..#G_2]];
    Append(~H_2, elt);
  end for;
  C_2222 := [J_p!(k*D2+l*D4+m*D5+n*D21) : k,l,m,n in [0,1]];
  for j in [1..#C_2222] do
    if C_2222[j] in H_2 then
      found[j,2] := true;
    end if;
  end for;
  print p,"Torsion <= 16 ?",&and[f[2] : f in found],Sort([Order(g) : g in Generators(G)]);
end for;


// --- rank possibly two (< genus so Chabauty might apply)
print "\n------Rank------";
//[[i,j,k,l,m,n] : i,j in [-4..4],k,l,m,n in [0,1],d in D|  d eq i*D0+j*D1+k*D2+l*D4+m*D5+n*D21];
print TwoSelmerGroup(J3);