////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////
//    Programs related to the                             //
//    computations of hyperelliptic logarithms            //
//    and Periods of                                      //
//    Hyperelliptic curves                                //
//    Started 28 June 2010                                //
//    Homero R. Gallegos-Ruiz                             //
//                                                        //
////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////

// Given six real numbers [a,a1,b,b1,c,c1] in ascending order
// and an integer value for precision
// this function computes to the requiered precision
// twice the integrals from a to a1,
// from b to b1, and from c to c1 of 1/Sqrt(-P(x))
// and of x/Sqrt(-P(x))
// where P(x)=(x-a)(x-a1)(x-b)(x-b1)(x-c)(x-c1);
// See Apendix A2 of Bost & Mestre
// The output is a 2x3 matrix, each column having
// the integrals corresponding to the pair of points 
// a a1, b b1, c c1;
// RF, real field where everything is defined
complexPeriods:=function(rts,precision,RF);
a:=RF!rts[1];
a1:=RF!rts[2];
b:=RF!rts[3];
b1:=RF!rts[4];
c:=RF!rts[5];
c1:=RF!rts[6];
D:=RF!1;
for i:=1 to precision do; // To the moment, I haven't computed the precise number the steps that will
//give the integral to the required precision, but the following works
//Each step 'doubles precision', more precisely, the convergence of the sequences ai,bi,ci is quadratic.
//Then the for has more than enough steps to get the desired precision
// The true integral int_a^a'(dx/y) is in between 2*Sqrt(Dn) * Pi / ((a-b1)*(a-c1)) and 2*Sqrt(Dn) * Pi / ((a1-b)*(a1-c))
// so one of them should give the integral to the required precision if the difference between the two is < E(-precision)
stepPrecision:=[2*Sqrt(D) * Pi(RF) *(1/ ((a1-b)*(a1-c)) - 1/((a-b1)*(a-c1))), 2*Sqrt(D) * Pi(RF) *(1/ ((b-a1)*(c-b1)) - 1/((b1-a)*(c1-b))), 2*Sqrt(D) * Pi(RF) *(1/ ((c-a1)*(c-b1)) - 1/((c1-a)*(c1-b))), 2*Sqrt(D) * Pi(RF) *(a1/ ((a1-b)*(a1-c)) - a/((a-b1)*(a-c1))), 2*Sqrt(D) * Pi(RF) *(b1/ ((b-a1)*(c-b1)) - b/((b1-a)*(c1-b))), 2*Sqrt(D) * Pi(RF) *(c1/ ((c-a1)*(c-b1)) - c/((c1-a)*(c1-b))) ];
stepPrecision:=Max(stepPrecision); //This variable is telling us the maximum of the difference between
//the bounds for every one of the six integrals.
//print stepPrecision;  //Print the precision at each step.
if stepPrecision lt 10^(-precision) then;
print "Precision attained in ", i, " steps."; // The number of steps will be printed
break i;  // This will stop the computations to give the integral to the required precision.
end if;
A:=(Sqrt((c - b)*(c - b1)*(c1 - b)*(c1 - b1)));
B:=(Sqrt((c - a)*(c - a1)*(c1 - a)*(c1 - a1)));
C:=(Sqrt((b - a)*(b - a1)*(b1 - a)*(b1 - a1)));
u:=((c*c1 - a*a1 - B)/(c + c1 - a - a1));
u1:=((b*b1 - a*a1 - C)/(b + b1 - a - a1));
v:= ((b*b1 - a*a1 + C)/(b + b1 - a - a1)); 
v1:= ((c*c1 - b*b1 - A)/(c + c1 - b - b1));
w:= ((c*c1 - b*b1 + A)/(c + c1 - b - b1)); 
w1:= ((c*c1 - a*a1 + B)/(c + c1 - a - a1));
D:=(4*D*(a*a1*(c + c1 - b - b1) + b*b1*(a + a1 - c - c1) + \
  c*c1*(b + b1 - a - a1))/((b + b1 - a - a1)*(c + c1 - b - b1)*(c + \
    c1 - a - a1)));
a := u; a1 := u1; b:= v; b1 := v1; c := w; c1 := w1;
end for;
V1:=(2*(Pi(RF)*Sqrt(D))/((a1 - b)*(a1 - c))); //This is the integral
			       // of 1/Sqrt(P(x)) from a to a1;
V2:=(2*(Pi(RF)*Sqrt(D))*a1/((a1 - b)*(a1 - c))); //This is the integral
			       // of x/Sqrt(P(x)) from a to a1
W1:=(2*(Pi(RF)*Sqrt(D))/((b - a1)*(c-b1))); // and so on.
W2:=(2*(Pi(RF)*Sqrt(D))*b1/((b - a1)*(c-b1)));
U1:=(2*(Pi(RF)*Sqrt(D))/((c - b1)*(c - a1)));
U2:=(2*(Pi(RF)*Sqrt(D))*c1/((c - b1)*(c - a1)));
N:=Matrix(2,3,[V1,W1,U1,V2,W2,U2]);
return N;
end function;

// This function computes the image of the roots of a polynomial
// under a linear transformation that changes their order
// INPUT an ordered list with six different real elements 
// (the real roots of a degree six polynomial)
//  defined over RF, a real field
invertedRts:=function(rts,RF);
a:=RF!rts[1];
ap:=RF!rts[2];
return [1/(2*rts[1]-a-ap),1/(2*rts[6]-a-ap),1/(2*rts[5]-a-ap),1/(2*rts[4]-a-ap),1/(2*rts[3]-a-ap),1/(2*rts[2]-a-ap)];
end function;


// This function computes the real integrals
realPeriods:=function(rts,precision,RF);
a:=RF!rts[1];
a1:=RF!rts[2];
rts:=invertedRts(rts,RF); //Following Bost&Mestre, we change the basis
c:=Sqrt(-&*rts);       //and compute the matrix of periods of the new curve
periods:=c*complexPeriods(rts,precision,RF);  
L:=Matrix(2,2,[RF!0,4,2,2*(a+a1)]); //This matrices will change the integrals
R:=Matrix(3,3,[RF!0,0,1,0,1,0,1,0,0]); //to the matrices wrt current coordinates
return L*periods*R;
end function;

// Instead of giving the roots, now we give a polynomial with rational coeffs
// Input: f a degree 6 polynomial MONIC, with rational coefficients
// with real roots only, all different
// precision: The desired precision for the integrals
periodsRealCase:=function(f, precision);
F:=RealField(3*precision); //We want to compute quotients, square roots and
			   //many other things, we need even a higher precision
SetDefaultRealField(F);    //This saves some work later, as we don't have to 
			   //redefine or coerce elements into fields of
			   // higher precision
x:=Parent(f).1;
Rx<x>:=PolynomialRing(RealField()); //We see now the polynomial as defined
f:=Rx!f; 	       	            //over the reals
rts:=[Roots(f)[1,1],Roots(f)[2,1],Roots(f)[3,1],Roots(f)[4,1],\
	   Roots(f)[5,1],Roots(f)[6,1]];
rts:=Sort(rts);
PeriodMatrix:=realPeriods(rts,precision,F);
return PeriodMatrix;
end function;

// We now want to compute integrals of the form 
// int_a^x, int_b^x, int_c^x with x in [a,a'], [b,b'] or [c,c']

// The following function computes a list of values x_1, dots x_{2^n}, 
// a corresponding list of signs e_i, a factor D, and three real numbers f1,f2,f3 such that 
// the desired integral equals
// f1 Ia + f2 Ib + f3 Ic + D * sum e_i int_?^{x_i} 
// where ? equals a_{n+1}, b_{n+1}, or c_{n+1} depending on the location of the x_i
// Input: d, the factor from the previous step
// fr: the three real numbers from the previous step
// lx: the list of values and signs from the previous step
// rts: a_n, a_n', b_n, b_n', c_n, c_n'
// n, an integer, the step number
// RF, the real field where everything is defined.
// needed as I don't grasp yet how MAGMA deals with real arithmetic
stepIntegral:=function(d, fr, lx, rts, n, RF);
a:=RF!rts[1];
a1:=RF!rts[2];
b:=RF!rts[3];
b1:=RF!rts[4];
c:=RF!rts[5];
c1:=RF!rts[6];
// a,a1,b,b1,c,c1:=Explode(rts);
fra, frb, frc:= Explode(fr);
newlx:=[];
newd:=(4*d*(a*a1*(c + c1 - b - b1) + b*b1*(a + a1 - c - c1) + \
  c*c1*(b + b1 - a - a1))/((b + b1 - a - a1)*(c + c1 - b - b1)*(c + \
  c1 - a - a1)));  // This is the T sub n in BM
A:=(Sqrt((c - b)*(c - b1)*(c1 - b)*(c1 - b1)));  // Discriminants of [Q,R],
B:=(Sqrt((c - a)*(c - a1)*(c1 - a)*(c1 - a1)));  // [R,P]
C:=(Sqrt((b - a)*(b - a1)*(b1 - a)*(b1 - a1)));  // and [P,Q]
newrts:=[((c*c1 - a*a1 - B)/(c + c1 - a - a1)),
((b*b1 - a*a1 - C)/(b + b1 - a - a1)),
((b*b1 - a*a1 + C)/(b + b1 - a - a1)),
((c*c1 - b*b1 - A)/(c + c1 - b - b1)),
((c*c1 - b*b1 + A)/(c + c1 - b - b1)),
((c*c1 - a*a1 + B)/(c + c1 - a - a1))]; // Base points in the next curve
for xi in lx do;
     if xi[1] le a1 then
       phi0:=a*a1*b + a*a1*b1 - a*a1*c - a*a1*c1 - a*b*b1 - a1*b*b1 + b*b1*c + b*b1*c1+
        xi[1]*(a*c + a*c1 + a1*c + a1*c1 - b*c - b*c1 - b1*c - b1*c1)+
        xi[1]^2*(-a - a1 + b + b1);
       phi1:= 2*a*a1*c*c1 - 2*b*b1*c*c1+
        xi[1]*(-2*a*a1*b - 2*a*a1*b1 + 2*a*b*b1 - 2*a*c*c1 + 2*a1*b*b1 - 2*a1*c*c1 + 2*b*c*c1 + 2*b1*c*c1)+
        xi[1]^2*(2*a*a1 - 2*b*b1);
       det:=4*(xi[1]-c)*(xi[1]-c1)*(xi[1]*(a*a1 - a*c1 - a1*c1 - b*b1 + b*c1 + b1*c1)
        - a*a1*b - a*a1*b1 + a*a1*c1 + a*b*b1 + a1*b*b1 - b*b1*c1)*
        (xi[1]*(a*a1 - a*c - a1*c - b*b1 + b*c + b1*c) - a*a1*b - a*a1*b1 + a*a1*c + a*b*b1 + a1*b*b1 - b*b1*c);
       z1:=(-phi1-Sqrt(det))/(2*phi0);
       z2:=(-phi1+Sqrt(det))/(2*phi0);
          if xi[1] le newrts[2] then            // The sign of y in the new curve changes at x=newrts[2]
              frc:=frc-xi[2]/2^n;               // See the paper to see why are we adding this.
              newlx:=Append(newlx,<RF!z1,xi[2]>);
              newlx:=Append(newlx,<RF!z2,xi[2]>);              
          else
              fra:=fra+xi[2]/2^(n-1);
              frc:=frc-xi[2]/2^n;
              newlx:=Append(newlx,<RF!z1,-xi[2]>);
              newlx:=Append(newlx,<RF!z2,xi[2]>); 
          end if;
     elif xi[1] le b1 then
        phi0:= -a*b*b1 + a*c*c1 - a1*b*b1 + a1*c*c1 + b*b1*c + b*b1*c1 - b*c*c1 - b1*c*c1+
            xi[1]*(a*b + a*b1 - a*c - a*c1 + a1*b + a1*b1 - a1*c - a1*c1)+
            xi[1]^2*(-b - b1 + c + c1);
        phi1:= 2*a*a1*b*b1 - 2*a*a1*c*c1+
            xi[1]*(-2*a*a1*b - 2*a*a1*b1 + 2*a*a1*c + 2*a*a1*c1 - 2*b*b1*c - 2*b*b1*c1 + 2*b*c*c1 + 2*b1*c*c1)+
            xi[1]^2*( 2*b*b1 - 2*c*c1);
        det:=4*(xi[1]-a)*(xi[1]-a1)*(xi[1]*(a1*b + a1*b1 - a1*c - a1*c1 - b*b1 + c*c1)
            - a1*b*b1 + a1*c*c1 + b*b1*c + b*b1*c1 - b*c*c1 - b1*c*c1)*
            (xi[1]*(a*b + a*b1 - a*c - a*c1 - b*b1 + c*c1)-a*b*b1 + a*c*c1 + b*b1*c + b*b1*c1 - b*c*c1 - b1*c*c1);
z1:=(-phi1+Sqrt(det))/(2*phi0);   // Here the leading coefficient phi0 is negative, hence the order of the roots changes.
        z2:=(-phi1-Sqrt(det))/(2*phi0);
          if xi[1] le newrts[4] then
              fra:=fra+xi[2]/2^n;
              newlx:=Append(newlx,<RF!z1,-xi[2]>);
              newlx:=Append(newlx,<RF!z2,xi[2]>);              
          else
              frb:=frb+xi[2]/2^(n-1);
              fra:=fra+xi[2]/2^n;
              newlx:=Append(newlx,<RF!z1,-xi[2]>);
              newlx:=Append(newlx,<RF!z2,-xi[2]>); 
          end if;
     else
        phi0:=  a*a1*b + a*a1*b1 - a*a1*c - a*a1*c1 + a*c*c1 + a1*c*c1 - b*c*c1 - b1*c*c1+
           xi[1]*(-a*b - a*b1 - a1*b - a1*b1 + b*c + b*c1 + b1*c + b1*c1)+
           xi[1]^2*(a + a1 - c - c1);
        phi1:= -2*a*a1*b*b1 + 2*b*b1*c*c1+
           xi[1]*(2*a*a1*c + 2*a*a1*c1 + 2*a*b*b1 - 2*a*c*c1 + 2*a1*b*b1 - 2*a1*c*c1 - 2*b*b1*c - 2*b*b1*c1)+
           xi[1]^2*(-2*a*a1 + 2*c*c1);
        det:=4*(xi[1]-b)*(xi[1]-b1)*(xi[1]*(a*a1 - a*b1 - a1*b1 + b1*c + b1*c1 - c*c1)
           + a*a1*b1 - a*a1*c - a*a1*c1 + a*c*c1 + a1*c*c1 - b1*c*c1)*
           (xi[1]*(a*a1 - a*b - a1*b + b*c + b*c1 - c*c1)+a*a1*b - a*a1*c - a*a1*c1 + a*c*c1 + a1*c*c1 - b*c*c1);
        z1:=(-phi1+Sqrt(det))/(2*phi0);
        z2:=(-phi1-Sqrt(det))/(2*phi0);
          if xi[1] le newrts[6] then
              frb:=frb+xi[2]/2^n;
              newlx:=Append(newlx,<RF!z1,-xi[2]>);
              newlx:=Append(newlx,<RF!z2,xi[2]>);              
          else
              frc:=frc+xi[2]/2^(n-1);
              frb:=frb+xi[2]/2^n;
              newlx:=Append(newlx,<RF!z1,-xi[2]>);
              newlx:=Append(newlx,<RF!z2,-xi[2]>); 
          end if;
     end if;
end for;
newfr:=[fra,frb,frc];
return newd, newfr, newlx, newrts;
end function;

// The following function computes an upper bound and a lower bound for the desired integral,
// and the difference between the two to get a bound on the precision of the integral
// Same input as before
precisionBound:=function(d, fr, lx, rts, n, RF);
pi:=Pi(RF);
a:=RF!rts[1];
a1:=RF!rts[2];
b:=RF!rts[3];
b1:=RF!rts[4];
c:=RF!rts[5];
c1:=RF!rts[6];
//a,a1,b,b1,c,c1:=Explode(rts);
A:=[a,a1]; B:=[b,b1]; C:=[c,c1];
// The following are lower and upper bounds for the half periods of the n-th curve.
boundsa:=[1/((b1-a)*(c1-a)), 1/((b-a1)*(c-a1))];
boundsb:=[1/((b1-a)*(c1-b)), 1/((b-a1)*(c-b1))];
boundsc:=[1/((c1-b)*(c1-a)), 1/((c-b1)*(c-a1))];
fra,frb,frc:=Explode(fr); // The integral is given by fra Ia + frb Ib + frc Ic + the other integrals
upperbound:=<RF!0.0,RF!0.0>;
lowerbound:=<RF!0.0,RF!0.0>;
// We now compute upper and lower bounds individually for every element of the final sum
// First, the portions of the half integral.
// We are computing simultaneously the integrals int dx/y and x*dx/y.
// We store the respective bounds in tuples <b1,b2> where b1 is the bound for the first
// differential and b2 the bound for the second.
if fra ne 0 then  //If fra eq 0 the following computation will fail. Hence the if.
   upperbound:= upperbound + < fra*pi*Sqrt(d)*boundsa[Integers()!(Sign(fra)/2+3/2)], \
   fra*pi*Sqrt(d)*A[Integers()!(Sign(fra)/2+3/2)]*boundsa[Integers()!(Sign(fra)/2+3/2)] >;  
                                                    // When the sign of fra is negative
                                                    // the lower bound for Ia gives an upper bound for 
                                                    // fra*Ia. The function x/2+3/2 takes care of the entry
                                                    // we need, 2 if x=1, 1 if x=-1.
   lowerbound:= lowerbound + <fra*pi*Sqrt(d)*boundsa[Integers()!(-Sign(fra)/2+3/2)], \
   fra*pi*Sqrt(d)*A[Integers()!(-Sign(fra)/2+3/2)]*boundsa[Integers()!(-Sign(fra)/2+3/2)]>;
end if;
if frb ne 0 then
   upperbound:= upperbound + < frb*pi*Sqrt(d)*boundsb[Integers()!(Sign(frb)/2+3/2)], \
   frb*pi*Sqrt(d)*B[Integers()!(Sign(frb)/2+3/2)]*boundsb[Integers()!(Sign(frb)/2+3/2)] >; 
   lowerbound:= lowerbound + <frb*pi*Sqrt(d)*boundsb[Integers()!(-Sign(frb)/2+3/2)], \
   frb*pi*Sqrt(d)*B[Integers()!(-Sign(frb)/2+3/2)]*boundsb[Integers()!(-Sign(frb)/2+3/2)]>;
end if;
if frc ne 0 then
   upperbound:= upperbound + < frc*pi*Sqrt(d)*boundsc[Integers()!(Sign(frc)/2+3/2)], \
   frc*pi*Sqrt(d)*C[Integers()!(Sign(frc)/2+3/2)]*boundsc[Integers()!(Sign(frc)/2+3/2)] >;  
   lowerbound:= lowerbound + <frc*pi*Sqrt(d)*boundsc[Integers()!(-Sign(frc)/2+3/2)], \
   frc*pi*Sqrt(d)*C[Integers()!(-Sign(frc)/2+3/2)]*boundsc[Integers()!(-Sign(frc)/2+3/2)]>;
end if;
// Now, individual bounds for every xi in the list, taking care of the signs.
for xi in lx do;
   if xi[1] le a1 then
       arc:=Arcsin(Sqrt((xi[1]-a)/(a1-a)));
       upperbound:=upperbound + < xi[2]*arc*Sqrt(d)/2^(n-2)*boundsa[Integers()!(xi[2]/2+3/2)] , \
          xi[2]*arc*Sqrt(d)/2^(n-2)*boundsa[Integers()!(xi[2]/2+3/2)]*A[Integers()!(xi[2]/2+3/2)] >;
       lowerbound:=lowerbound + < xi[2]*arc*Sqrt(d)/2^(n-2)*boundsa[Integers()!(-xi[2]/2+3/2)] , \
          xi[2]*arc*Sqrt(d)/2^(n-2)*boundsa[Integers()!(-xi[2]/2+3/2)]*A[Integers()!(-xi[2]/2+3/2)] >;
   elif xi[1] le b1 then
       arc:=Arcsin(Sqrt((xi[1]-b)/(b1-b)));
       upperbound:=upperbound + < xi[2]*arc*Sqrt(d)/2^(n-2)*boundsb[Integers()!(xi[2]/2+3/2)] , \
          xi[2]*arc*Sqrt(d)/2^(n-2)*boundsb[Integers()!(xi[2]/2+3/2)]*B[Integers()!(xi[2]/2+3/2)] >;
       lowerbound:=lowerbound + < xi[2]*arc*Sqrt(d)/2^(n-2)*boundsb[Integers()!(-xi[2]/2+3/2)] , \
          xi[2]*arc*Sqrt(d)/2^(n-2)*boundsb[Integers()!(-xi[2]/2+3/2)]*B[Integers()!(-xi[2]/2+3/2)] >;
   else
       arc:=Arcsin(Sqrt((xi[1]-c)/(c1-c)));
       upperbound:=upperbound + < xi[2]*arc*Sqrt(d)/2^(n-2)*boundsc[Integers()!(xi[2]/2+3/2)] , \
          xi[2]*arc*Sqrt(d)/2^(n-2)*boundsc[Integers()!(xi[2]/2+3/2)]*C[Integers()!(xi[2]/2+3/2)] >;
       lowerbound:=lowerbound + < xi[2]*arc*Sqrt(d)/2^(n-2)*boundsc[Integers()!(-xi[2]/2+3/2)] , \
          xi[2]*arc*Sqrt(d)/2^(n-2)*boundsc[Integers()!(-xi[2]/2+3/2)]*C[Integers()!(-xi[2]/2+3/2)] >;
   end if;
end for;
stepPrecision:=Max([upperbound[1]-lowerbound[1],upperbound[2]-lowerbound[2]]);
return stepPrecision, lowerbound;
end function;

// Now we should be ready for a simple version of the computation of 
// int from a to x of dx/y, xdx/y with x in [a,a']
// and the corresponding version for [b,b'], [c,c']
// input: rts, six different real numbers at very high precision, 
// x1 a real value lying on either [a,a'], [b,b'], [c,c']
// this condition is not checked, we only assume it is fulfilled.
// precision, the exponent of the desired precision of computation for the integrals
// RF, the real field where everything is defined
// needed, as I don't grasp how MAGMA deals with real arithmetic
simplifiedHyperLog:=function(rts, x1, precision, RF);
lx:=[<RF!x1,1>];  
d:=RF!1;   // A factor relating the integrals
fr:=[0,0,0];  //  Portions of the half periods we will need
for i:=1 to precision do;  // We don't need many steps, but I don't have a better option to check the precision
stepPrecision,lowerbound:=precisionBound(d,fr,lx,rts,i,RF);    // see the description of precisionBound
   if stepPrecision lt 10^(-precision) then  // here we check the precision
   print "Precision required attained in ", i, " steps.";  // Only for information purposes.
      break i;
   end if;
d, fr, lx, rts:=stepIntegral(d,fr,lx,rts,i,RF);  // This is the integral after the change of variables
end for;
// stepPrecision;  // Only for information purposes at the time of building the function
return lowerbound;
end function;

// Let f be a degree 5 polynomial defined over the rationalas,
// with 5 real roots, monic.
// Let e5 be the largest of its roots.
// This function computes twice the integral from e5 to infinity
// of dx/Sqrt(f(x)) and of xdx/Sqrt(f(x)).
// input: f, a degree 5 polynomial with rational coefficients,
// we assume it has 5 different real roots
// precision: an integer, the desired precision
PeriodsDegreeFive:=function(f,precision);
F:=RealField(3*precision); //We want to compute quotients, square roots and
			   //many other things, we need even a higher precision
SetDefaultRealField(F);    //This saves some work later, as we don't have to 
			   //redefine or coerce elements into fields of
			   // higher precision
x:=Parent(f).1;
Rx<x>:=PolynomialRing(RealField()); //We see now the polynomial as defined
f:=Rx!f; 	       	            //over the reals
rts:=[Roots(f)[1,1],Roots(f)[2,1],Roots(f)[3,1],Roots(f)[4,1],\
	   Roots(f)[5,1]];
rts:=Sort(rts);
e1,e2,e3,e4,e5:=Explode(rts);
newrts:=[1/(e2-(e2+e3)/2),1/(e1-(e2+e3)/2),1/(e5-(e2+e3)/2),\
            1/(e4-(e2+e3)/2),1/(e3-(e2+e3)/2)];   // These are the nonzero
           //  roots of the degree 6 polynomial that results from the
           //  change of variables mentioned in the text.
constantterm:=Sqrt(&*newrts); // This factor comes from the change of variables
newrts:=Insert(newrts,3,0);   // I didn't want to put 0 in the
                              // previous product! But 0 is the sixth
                              // root of the new polynomial
MatrixOfPeriods:=complexPeriods(newrts,precision, F);  // Bost &  Mestre
  // algorithm for the computation of the periods with respect to the degree
  // 6 curve.
  // The integral e1 to e2 is encoded in the first column, with the
  //  oposite signs, the integral e3 to e4 is encoded in the last
  //  column, and the integral e5 to infinty is encoded in the middle
  //  column. As the first column minus the second plus the third is
  //  0, we will only keep the first 2 columns.
period1:=constantterm*MatrixOfPeriods[2];   // These are the integrals
                                            // of 1/Sqrt(f(x))
                        // with respect to the degree 5 curve. 
period2:=constantterm*(MatrixOfPeriods[1]+(e2+e3)/2*MatrixOfPeriods[2]);
                     // These are the integrals of x/Sqrt(f(x))
p1:=[-period1[1], -period2[1]];
p2:=[period1[2], period2[2]];
// We need to change the signs of the first column as noted before.
// We store the first column in p1, and the second in p2.  
return p1,p2;
end function;

// Let f be a degree 5 monic polynomial, defined over the rationals
// and with 5 real roots, e1 < e2 < e3 < e4 < e5.
// Let x0 be a real number lying either in [e1, e2], [e3, e4] or [e5, infty)
// This function computes the integral from x0 to ei, where ei is either e2,
// e4 or infinity, depending on the location of x1.
// Input: x0, a real number
// f, degree 5 polynomial monic with 5 real roots
// precision: an integer, the required digits of precision.
HyperellipticLogDegree5:=function(x0,f,precision);
F:=RealField(10*precision); //We want to compute quotients, square roots and
			   //many other things, we need even a higher precision.
SetDefaultRealField(F);    //This saves some work later, as we don't have to 
			   //redefine or coerce elements into fields of
			   // higher precision
x:=Parent(f).1;
Rx<x>:=PolynomialRing(RealField()); //We see now the polynomial as defined
f:=Rx!f; 	       	            //over the reals
rts:=[Roots(f)[1,1],Roots(f)[2,1],Roots(f)[3,1],Roots(f)[4,1],\
	   Roots(f)[5,1]];
rts:=Sort(rts);
e1,e2,e3,e4,e5:=Explode(rts);
newrts:=[1/(e2-(e2+e3)/2),1/(e1-(e2+e3)/2),1/(e5-(e2+e3)/2),\
            1/(e4-(e2+e3)/2),1/(e3-(e2+e3)/2)];   // These are the nonzero
           //  roots of the degree 6 polynomial that results from the
           //  change of variables mentioned in the text.
constantterm:=Sqrt(&*newrts); // This factor comes from the change of variables
newrts:=Insert(newrts,3,0);   // I didn't want to put 0 in the
                              // previous product! But 0 is the sixth
                              // root of the new polynomial
x1:=1/(x0-(e2+e3)/2);  // This is the value of x after the change of coordinates
logs:=simplifiedHyperLog(newrts, x1, precision, F);   // This are the
         // integrals from the nearest root in newrts that is less
	 // than x1.
log1:=constantterm*logs[2];   // These are the integrals
        // with respect to the degree 5 curve. 
log2:=constantterm*(logs[1]+ (e2+e3)/2*logs[2]);
if x0 le (e2+e3)/2 then;  // The value of the logs contains a AbsValue(X)
log1:=-log1; // the values le than (e2+e3)/2 are mapped to negative values
log2:=-log2; // hence the change of sign.
end if;
return log1,log2;
end function;
