//////////////////////////////////////////////////////////////
// Computations for the proof of Theorem 1.1 in             //
// Computing Integral Points on Genus 2 Hyperelliptic       //
// Curves Estimating Hyperelliptic Logarithms               //
// H.R. Gallegos-Ruiz                                       //
// February 2017                                            //
//////////////////////////////////////////////////////////////
Q<x>:=PolynomialRing(Rationals());
f:=x^5 + 105*x^4 + 4405*x^3 + 92295*x^2 + 965794*x + 4038280;
a5:=LeadingCoefficient(f);
C:=HyperellipticCurve(f);
print "Computing integral points of C: y2 = ", f;
J:=Jacobian(C);
r:=RankBounds(J);
print "The rank of the MW group for this curve is ", r;
T:=TorsionSubgroup(J);
print "The torsion subgroup is ", T;
t:=Exponent(T);
bas := ReducedBasis(Points(J : Bound :=1000));
a:=1;
binarycoefficients:=[];
print "Computing a full set of representatives for J/2J...";
// What we do here is finding all expressions of five 0s and 1s to
// multiply by the elements in bas and find the representatives for J/2J
for i in [0..31] do binarycoefficients:=Append(binarycoefficients,
     Prune(Intseq(i+32,2))); end for;
fullSetReps:=[];
for i in [1..32] do 
  fullSetReps:=Append(fullSetReps,binarycoefficients[i][1]*bas[1]
    + binarycoefficients[i][2]*bas[2] + binarycoefficients[i][3]*bas[3]
    + binarycoefficients[i][4]*bas[4] + binarycoefficients[i][5]*bas[5]); 
end for;
print "done.";
print "Computing Kappa set...";
// The following computations are in the context of Prop 7.1 and theorem 8.1
// of my paper "S-integral points on hyperelliptic curves".
K<alpha>:=NumberField(f);
kappaset:=[];
for i in fullSetReps do
      kappa:=Evaluate(i[1], alpha);  
      d1,d2:=SquarefreeFactorization(Denominator(kappa));
      d:=d1*d2;
      kappa:=d^2*i[1];
      if i[3] eq 1 then kappa:=a*kappa; end if;
      kappaset:=Append(kappaset,kappa);
end for;
print "done.";
load "upperbounds.m";
S:=[];  // We are dealing with integral points.
Bounds:=[];
print "Computing an upper bound for the integral points on C...";
time for k in kappaset do;
  Bounds:= Append(Bounds,heightxBoundx(f,k,S));
end for;
upperbound:= Max(Bounds);
print "done.";
print "An integral point on C has log x at most ", upperbound;
// Compute a lower bound for the height difference
load "heightdifference.m";
mu1:=HeightDifferenceLowerBound(f);
print "A lower bound for the high differente h - hhat is ", mu1;
mu2:= Min([l[1] : l in Eigenvalues(HeightPairingMatrix(bas))]);
print "The smallest eigenvalue of the height pairing matrix is ", mu2;
// Compute a bound for the coefficients in  n1 D1+ ... + n5 D5
// Corollary 3.2
N:= Sqrt(mu2^(-1)*(Log(a5) + 2*upperbound-mu1));
print "an upper bound for the norm of (n1,...,n5) is ", N;
//  The choices for K are based in the following quantity
//  ((48/5)*(r*N*t + (1/4)*Sqrt(r)*N + (5/48)*N + 1))^((r+4)/4);
print "Setting up reduction process...";
K:=10^970;
F:=RealField(1300);
SetDefaultRealField(F);
print "Computing big period matrix...";
time AJ:=AnalyticJacobian(C);
M:=BigPeriodMatrix(AJ);
RealM:=Matrix(4, [Re(M[1,1]), Re(M[1,2]), Re(M[1,3]), Re(M[1,4]), Im(M[1,1]), Im(M[1,2]), Im(M[1,3]), Im(M[1,4]), Re(M[2,1]), Re(M[2,2]), Re(M[2,3]), Re(M[2,4]),  Im(M[2,1]), Im(M[2,2]), Im(M[2,3]), Im(M[2,4])]);
print "done.";
R<x>:=PolynomialRing(ComplexField());
print "Computing hyperelliptic logarithms of the generators of the MW group...";
analyticPts:=[];
// The roots of the first polynomial on the divisor representation are
// the x coordinates for the points on the support, the y coordinates
// are obtained in terms of x by the second polynomial
for pt in bas do;
   coords := [r[1] : r in Roots(Evaluate(pt[1],x))];
   coords := [<d1,Evaluate(pt[2],d1)> : d1 in coords];
analyticPts:=Append(analyticPts,coords);
end for;
hyperlogs:=[];
// ToAnalyticJacobian computes the hyperelliptic logarithm modulo Lambda
for pt in analyticPts do;
   hyperlogpt:= &+[ToAnalyticJacobian(d[1],d[2],AJ): d in pt];
   hyperlogpt:=[Re(hyperlogpt[1][1]), Im(hyperlogpt[1][1]), Re(hyperlogpt[2][1]), Im(hyperlogpt[2][1])];
   coordsInOmegas:=Eltseq(RealM^(-1)*Matrix(1,hyperlogpt));
   coordsReduced:=[];  // We want to find the hyper log contained in the fundamental domain
   for i in coordsInOmegas do;
     coordsReduced:=Append(coordsReduced, i-Floor(i));
   end for;
   hyperlogpt:=RealM*Matrix(1,coordsReduced);
   hyperlogs:=Append(hyperlogs,hyperlogpt);
end for;
print "done.";
// Block matrices to set up matrix A depending on K.
B1:=ScalarMatrix(r,1);
B2:=ZeroMatrix(Integers(),r,4);
B3:=HorizontalJoin(B1,B2);
B4:=K*Transpose(Matrix(hyperlogs));
B5:=K*RealM/t;
B6:=HorizontalJoin(B4,B5);
D:=Eltseq(B6);
B7:=[];
for i in D do; B7:=Append(B7,Round(i)); end for;
B7:=Matrix(r+4,B7);
A:=VerticalJoin(B3,B7);
// MAGMA takes rows as vectors generating the lattice, in the paper we take columns.
L:=Lattice(Transpose(A));
l:=Sqrt(Minimum(L));
assert(l gt ((1/5)*(48*Sqrt(r)*N*t + 12*Sqrt(r)*N + 5*N + 48)));  // When the condition is untrue, the assertion fails and the program stops
l;
((1/5)*(48*Sqrt(r)*N*t + 12*Sqrt(r)*N + 5*N + 48));
// Constants in Corollary 5.2:
K1:=2^(7/2)*Exp(-mu1/4)/Sqrt(a5);
K2:=mu2/4;
// Reduced bound according to Proposition 6.2
N:=Sqrt((Log(4*K*K1)- Log(l-((1/5)*(48*Sqrt(r)*N*t + 12*Sqrt(r)*N + 5*N + 48))))/K2);
print "the new bound for the coefficients wrt MW basis is ", N;
K:=10^13;
B4:=K*Transpose(Matrix(hyperlogs));
B5:=K*RealM/t;
B6:=HorizontalJoin(B4,B5);
D:=Eltseq(B6);
B7:=[];
for i in D do; B7:=Append(B7,Round(i)); end for;
B7:=Matrix(r+4,B7);
A:=VerticalJoin(B3,B7);
L:=Lattice(Transpose(A));
l:=Sqrt(Minimum(L));
assert(l gt ((1/5)*(48*Sqrt(r)*N*t + 12*Sqrt(r)*N + 5*N + 48)));
N:=Sqrt((Log(4*K*K1)- Log(l-((1/5)*(48*Sqrt(r)*N*t + 12*Sqrt(r)*N + 5*N + 48))))/K2);
print "the new bound for the coefficients wrt MW basis is ", N;
K:=2*10^10;
B4:=K*Transpose(Matrix(hyperlogs));
B5:=K*RealM/t;
B6:=HorizontalJoin(B4,B5);
D:=Eltseq(B6);
B7:=[];
for i in D do; B7:=Append(B7,Round(i)); end for;
B7:=Matrix(r+4,B7);
A:=VerticalJoin(B3,B7);
L:=Lattice(Transpose(A));
l:=Sqrt(Minimum(L));
assert(l gt ((1/5)*(48*Sqrt(r)*N*t + 12*Sqrt(r)*N + 5*N + 48)));
N:=Sqrt((Log(4*K*K1)- Log(l-((1/5)*(48*Sqrt(r)*N*t + 12*Sqrt(r)*N + 5*N + 48))))/K2);
print "the new bound for the coefficients wrt MW basis is ", N;
N0:=Floor(N);
pts:=[];
print "Computing all points of the form n1 D1+ ... + n5 D5";
print "with Euclidean norm at most N"; // Takes a long time
// We only compute those with i1 ge 0, as the others come from negatives.
time for i1:=0 to N0 do;
  for i2:=-N0 to N0 do;
    for i3:=-N0 to N0 do;
      for i4:=-N0 to N0 do;
        for i5:=-N0 to N0 do;
           if (i1^2 + i2^2 + i3^2 + i4^2 + i5^2 le N^2) then;
              pp:=i1*bas[1]+i2*bas[2]+i3*bas[3]+i4*bas[4]+i5*bas[5];
              if pp[3] eq 1 then;  // that means we have only one point P,
                pts:=Append(pts,pp);
              end if;
           end if;
         end for;
       end for;
     end for;
  end for;
end for;
print "The rational points found in C below the bound";
print "(up to negatives) are ", pts;
// Compute points with x not in the region covered by Proposition 6.2
print "Computing remaining points...";
boundx:=2*Max([Modulus(i[1]) : i in Roots(Evaluate(f,x))]);
boundx:=Max([boundx,1]);
ptsBelowBound:=Points(C: Bound:=Ceiling(boundx));
print "The rational points on C below the bound are ", ptsBelowBound;

