//////////////////////////////////////////////////////////////
// MAGMA script for the computations in                     //
// ``S-integral Points on Hyperelliptic Curves''            //
// by H. R. Gallegos-Ruiz                                   //
//////////////////////////////////////////////////////////////

// The references are numbered as in my paper
// ``S-integral Points on Hyperelliptic Curves''
// to appear in the International Journal of
// Number Theory.

// Consider the curve Y^2 - Y = X^5 - X and let S be the set of the
// first 22 primes.
// We will find the S-integral points of the curve. 
print "Computing the S-integral points on  Y^2 - Y = X^5 - X, ";
print "where S is the set of the first 22 rational primes.";

// First we find an upper bound for the size of the S-integral points
// on the curve 

load "upperbounds.m"; // Code for the computation of the upper bounds
                      // of the size of the integral points.

// Set up the curve.
// We transform the equation in to 
// C : 2*y^2 = x^5 - 16*x + 8
Qx<x>:=PolynomialRing(Rationals());
p:=x^5-16*x+8; // Polynomial after change of coordinates
               // y = 4*Y - 2; x = 2*X;

// The set K constructed in Lemma 2.1 is the following. 
kappaset:=[1,-2*x,4-2*x,-4-2*x,-2*x+x^2,2*x+x^2,-4+x^2,8*x-2*x^3]; 
                            // computed in [11,Proof of Theorem 1.1]
   S:=PrimesUpTo(NthPrime(22));
   print "S = primes up to" , NthPrime(22);
   Bounds:=[];
   time for k in kappaset do;
   	    print k;
           Bounds:= Append(Bounds,heightxBoundx(p,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;
   print "=======================";

// We now implement the variant of the Mordell--Weil sieve
// following [11] Yann Bugeaud, Maurice Mignotte, Samir Siksek, Michael Stoll
// and Szabolcs Tengely, "Integral points on hyperelliptic curves", Algebra and
// Number Theory 2 (2008) , no.8, 859--885 

load "MW.m"; // Code for the Mordell--Weil sieve

C:=HyperellipticCurve(4*x^5-4*x+1); // The curve after the change
                                    // of coordinates x = X; y = 2Y - 1
J:=Jacobian(C);
inf:=PointsAtInfinity(C)[1];
mwBasis:=[J![x,1],J![x-1,1],J![x^2-1,1]];    //Mordell--Weil basis for J

PC:=Points(C : Bound:=100); // known rational points
W:=[P-inf : P in PC];  // their image on the Jacobian

// From [11] we know that j(C) is contained in W + B*J, where
B:=4449329780614748206472972686179940652515754483274306796568214048000;

F:=FreeAbelianGroup(3);
L:=sub<F | [B*F.1,B*F.2,B*F.3]>;

// In terms of the map phi,  j(C(QQ)) is contained in W+phi(L);

// We now compute the constants in Lemma 9.1
mu1:=HeightConstant(J: Factor:=true, Effort:=2); mu1:=-mu1;
    // See '? HeightConstant' for the change of sign.
mu2:=Max([Sqrt(Height(i)):i in W]);
mu3:= Min([Sqrt(l[1]) : l in Eigenvalues(HeightPairingMatrix(mwBasis))]);
print "mu1 is ", mu1;
print "mu2 is ", mu2;
print "mu3 is ", mu3;

// we apply [8, Lemma 11.1] successively
// with the first 50000 primes (if of good reduction)
// Then, we improve B and start again.
// Four times are enough for contradicting the upper bounds
// for the set of $S$-integral points.

bprimes:=BadPrimes(J);
iprime:=2;  //initial prime
    for i:=1 to 4 do;
        time L,B,iL,mL,successLevels,bprimes,lprime:=mwSievePrimesEnhance(F,L,W,J,mwBasis,iprime,5*10^4,B,bprimes);
        print "Now j(C) is contained in W + newB times J where newB is";
        print B;
        assert(mu3*mL-mu2 gt 0); //Need to check according to lemma 9.1
        lowerBound:=(mu3*mL-mu2)^2+mu1;
        print "A lower bound for the size of j(P)";
        print "for P in the set of unknown rational points is";
        lowerBound;
    end for;

lbound:=(lowerBound-Log(2))/2;
print "If P=(x,y) is not one of the known rational points of 2y^2 = x^5 - 16x +8";
print "then the height of x is at least", lbound;
assert(lbound gt upperbound);
