// Computation of the integral points of the curve
// C: y^2 = x^5-5*x^3-1*x^2+3*x+1

// We first find an upper bound for the size of the integral points on C
// as in ``S-integral points on hyperelliptic curves''.
load "upperbounds.m";
Qx<x>:=PolynomialRing(Rationals());
f:=x^5-5*x^3-1*x^2+3*x+1;  // Our polynomial
C:=HyperellipticCurve(f);  
J:=Jacobian(C);
bas:=ReducedBasis(Points(J:Bound:=20));   // generators for the Mordell-Weil group
reps:=[J!0,bas[1],bas[2],bas[1]+bas[2]]; // representatives for J(Q)/2J(Q)
kappaset:=[];
for i in reps do;
  kappaset:=Include(kappaset, (-1)^(i[3])*i[1] );
end for;
S:=[]; // S is a finite set of primes when we compute S-integral
       // points. Here we need it to be empty.
Bounds:=[];                           //A place to save the bounds
time for k in kappaset do;
       Bounds:= Append(Bounds,heightxBoundx(f,k,S));  //Bound for each kappa
     end for;
upperbound:= Max(Bounds);
print "If P=(x,y) is an S-integral point of C, then";
print "h(x) is at most", upperbound;

Qy<y>:=PolynomialRing(RealField());  // I need to find 2 Max [AbsValue(alpha)], 
realf:=Evaluate(f,y);                // alpha running over the roots of f,
                                     // f has 5 real roots in this case
boundx:=2*Max([AbsoluteValue(i[1]) : i in Roots(realf) ]);

// Now we do the reduction method:
// We need first to setup some constants
// mu1, a lower bound for the Height Difference,
load "heightdifference.m";  // Code for the computation of the lower bound.
mu1:=HeightDifferenceLowerBound(f);

// mu3, the smallest eigenvalue of the Height Pairing Matrix
mu3:= Min([l[1] : l in Eigenvalues(HeightPairingMatrix(bas))]);

// K1, K2 real constants, such that the linear forms in 
// hyperelliptic logs is bounded by
// K1 Exp (-K2 M^2), where M is the bound for the coefficients
// of the linear form
K1:=(2^(7/2))*Exp(-mu1/4); // The leading coefficient of f is 1
K2:=mu3/4;

// K3, an upper bound for the coefficients of the linear form,
// We obtain it from the bound for the integral points
// computed before.
K3:=Sqrt(mu3^(-1)*(2*upperbound-mu1));
print "first upper bound for the coefficients of the linear form";
K3;

// Constants of Lemma 5.4.1
S:= 2*K3^2;
T:=(2*K3*(2+1) +3*2+1)/2;

// lower bound for the minimum of an LLL basis
// This is Theorem V.9 in [Smart]
// input M a square matrix with integer coefficients,
// This is assumed to be non singular.
// Output, a lower bound for the minimum of the lattice
// generated by the row vectors
// GramBound:=function(M);
// n:=Nrows(M);
// M:=MatrixAlgebra(RealField(),n)!M;
// vectors:=[M[j]: j in [1..n]];
// newvectors:=vectors;
// for i:=2 to n do;
//   for j:=1 to (i-1) do;
//     newvectors[i]:=newvectors[i] - (vectors[i],newvectors[j])/(newvectors[j],newvectors[j])*newvectors[j];
//   end for;
// end for;
// boundmin:=Max([Norm(vectors[1])/Norm(j) : j in newvectors   ]);
// return Sqrt(Norm(vectors[1])/boundmin);
// end function;

// We need to choose C about the size of K3^4
// we choose
C:=10^1000; 

load "magmaperiods1.m"; // Code for the computation of periods and
			// logs at required precision

p1,p2:=PeriodsDegreeFive(f,1000); // This are the vectors corresponding to
// the fundamental real periods.
l11,l12:=HyperellipticLogDegree5(-1.0,f,1000);  // This are half the hyperelliptic
l21,l22:=HyperellipticLogDegree5(0.0,f,1000);   // logs of 2Di with at least 1000 digits of precision
n11:=Round(C*l11);
n12:=Round(C*l21);
n13:=Round(C*p1[1]/2);
n14:=Round(C*p2[1]/2);
n21:=Round(C*l12);
n22:=Round(C*l22);
n23:=Round(C*p1[2]/2);
n24:=Round(C*p2[2]/2);
M:=Matrix(4,[1,0,n11,n21, 0,1,n12,n22, 0,0, n13,n23, 0,0, n14,n24]);
// LM:=LLL(M);
c4sq:=Min(Lattice(M));
assert(c4sq gt 2*T^2 + S);
K3:=Sqrt((Log(C*K1)-Log(Sqrt(c4sq-S)-T))/K2);
print "second upper bound for the coefficients of the linear form";
Floor(K3);

// This new upper bound is less than 90
// We then choose C to be
C:=10^8;
// Constants of Lemma 5.4.1
S:= 2*K3^2;
T:=(2*K3*(2+1) +3*2+1)/2;
n11:=Round(C*l11);
n12:=Round(C*l21);
n13:=Round(C*p1[1]/2);
n14:=Round(C*p2[1]/2);
n21:=Round(C*l12);
n22:=Round(C*l22);
n23:=Round(C*p1[2]/2);
n24:=Round(C*p2[2]/2);
M:=Matrix(4,[1,0,n11,n21, 0,1,n12,n22, 0,0, n13,n23, 0,0, n14,n24]);
//LM:=LLL(M);
c4sq:=Min(Lattice(M));
assert(c4sq gt 2*T^2 + S);
K3:=Sqrt((Log(C*K1)-Log(Sqrt(c4sq-S)-T))/K2);
print "third upper bound for the coefficients of the linear form";
Floor(K3);

// This new upper bound is less than 10
// We then choose C to be
C:=10^4;
// Constants of Lemma 5.4.1
S:= 2*K3^2;
T:=(2*K3*(2+1) +3*2+1)/2;
n11:=Round(C*l11);
n12:=Round(C*l21);
n13:=Round(C*p1[1]/2);
n14:=Round(C*p2[1]/2);
n21:=Round(C*l12);
n22:=Round(C*l22);
n23:=Round(C*p1[2]/2);
n24:=Round(C*p2[2]/2);
M:=Matrix(4,[1,0,n11,n21, 0,1,n12,n22, 0,0, n13,n23, 0,0, n14,n24]);
//LM:=LLL(M);
c4sq:=Min(Lattice(M));
assert(c4sq gt 2*T^2 + S);
K3:=Sqrt((Log(C*K1)-Log(Sqrt(c4sq-S)-T))/K2);
print "fourth upper bound for the coefficients of the linear form";
Floor(K3);
// This should be 7.  That is small enough


K3:=Floor(K3);
pts:=[];  // We will store in this list the divisors comming from the 
// linear combinations of the divisors up to the bound K3 that are of
// the form P-infty.  These are encoded in MAGMA as divisors of the 
// form (l(x), p(x),1)  The 1 means a degree 1 divisor, that is, one 
// point.  l(x) is a linear polynomial, its zero is the value of x we
// are looking for.  The corresponding point on C is
// (x, p(x))
for i:=-K3 to K3 do;
   for j:=-K3 to K3 do;
     pp:=i*bas[1]+j*bas[2];
     if pp[3] eq 1 then;  // that means we have only one point P,
        pts:=Append(pts,pp);
     end if;
   end for;
end for;
print "rational points below the bound", K3, "are";
pts;

//  We only need to look for points with 'small' x coordinate,
// Recall that the linear form is bounded as we claimed only
// for x > 2 Max(Norm(alphai)), alphai the roots of f
C:=HyperellipticCurve(f);  // We had changed the value of C
ptsC:=Points(C: Bound:=Ceiling(boundx));
print ptsC;
