/***
    This PARI/GP script is part of the following paper:

      P. J. Bruin, Bornes optimales pour la diff\'erence entre la
      hauteur de Weil et la hauteur de N\'eron--Tate sur les courbes
      elliptiques sur $\Qbar$.

    It contains an algorithm for computing the infimum and supremum of
    the local height difference function at the Archimedean place of
    an elliptic curve over \Q.

    The functions intended to be called directly are:

      phi_maximum(E, prec)
      phi_minimum(E, prec)
      phi_extrema(E, prec) = [phi_maximum(E, prec),
                              phi_minimum(E, prec)]

    The first two functions return a pair [z, a] with  z in \C
    and  a = \phi(z) such that

      \sup\phi - prec < a <= \sup\phi,  resp.
      \inf\phi <= a < \inf\phi + prec,

    where the \sup and \inf are taken over the period torus of  E.

    Example (output truncated):

      gp > \r hdiff.gp
      gp > E = ellinit("11a3");
      gp > phi_maximum(E, 0.001)
      %1 = [-4.536... - 0.005...*I, 0.5979737...]

    Thanks to Jan Steffen Müller for reporting several bugs.
***/

/*
  Local height function.  Alternative definition:

  Q(Ee, z) = real(Ee.C * z^2) + Ee.D * abs(z)^2

  lambda(E, Ee, z) = -real(ellsigma(E, z, 1)) + 1/2 * Q(Ee, z)
                     - 1/12 * log(abs(E.disc))
*/
lambda(om, z) = {
  local(tau = om[1]/om[2],
        w = z/om[2]);
  return(-log(abs(theta(exp(Pi*I*tau), Pi*w)/eta(tau, 1)))
         + Pi*imag(w)^2/imag(tau));
}

/*
  Return the two points in  \C/\Lambda  with given x-coordinate.
*/
ellxtoz(E, x) = {
  local(Y = ellordinate(E, x + 0.0*I),
        z = apply(y -> ellpointtoz(E, [x, y]), Y));
  return(if(length(z) == 1, [z[1], z[1]], z));
}

/*
  Two invariants of the period lattice of  E.
*/
lattice_C(E) = {
  return((E.eta[1]*conj(E.omega[2]) - E.eta[2]*conj(E.omega[1]))/(2*I*E.area));
}

lattice_D(E) = {
  return(Pi/E.area);
}

/*
  Initalise additional data for the elliptic curve  E.
*/
ellinit_extra(E) = {
  local(t = ellxtoz(E, 0),
        w = ellxtoz(E, 1)[1],
        int = -2*lambda(E.omega, w) + lambda(E.omega, w - t[1])
             + lambda(E.omega, w - t[2]),
        Ee = [lattice_C(E), lattice_D(E), t, int],
        M1, M2, K, W_bound);
  M1 = abs(lattice_C(E) + E.b2/12) + 1;
  M2 = abs(lattice_C(E) + E.b2/12) + (abs(E.b4) + abs(E.b6))/2;
  K = intnum(s = 0, 2*Pi,
             iferr(abs(4*exp(3*I*s) + E.b2*exp(2*I*s)
                       + 2*E.b4*exp(I*s) + E.b6)^(-1/2),
                   err, 0, errname(err) == "e_DOMAIN"));
  W_bound = max(abs(Zeta(E, Ee, w)) + M1 * K,
                abs(Zeta(E, Ee, w - t[1]) + Zeta(E, Ee, w - t[2]))/2 + M2 * K);
  Ee = concat(Ee, [max(M1, M2), W_bound]);
  return(Ee);
}

ee.C = ee[1];
ee.D = ee[2];
ee.t = ee[3];
ee.int = ee[4];
ee.M = ee[5];
ee.W_bound = ee[6];  \\ global bound on the function |W(z)|

/*
  A suitable modification of Weierstrass's zeta function,
  which is periodic but not meromorphic (because of the conj(z)).
*/
Zeta(E, Ee, z) = {
  return(ellzeta(E, z) - Ee.C * z - Ee.D * conj(z));
}

/*
  Check if the absolute value of the x-coordinate of the
  given point  z  in the torus corresponding to  E  is <= 1.
*/
xabs_le1(E, z) = {
  iferr(abs(ellwp(E.omega, z) - E.b2/12) <= 1, err, 0,
	errname(err) == "e_DOMAIN");
}

/*
  phi(z) is the local height difference function.
*/
phi(E, Ee, z) = {
  return(if(xabs_le1(E, z),
            -2*lambda(E.omega, z),
            -lambda(E.omega, z - Ee.t[1])
            - lambda(E.omega, z - Ee.t[2]) + Ee.int));
}

/*
  W(z) is the function such that
  d\phi(z) = W(z) dz + \bar W(z) d\bar z.
*/
W(E, Ee, z) = {
  return(if(xabs_le1(E, z),
            Zeta(E, Ee, z),
            (Zeta(E, Ee, z - Ee.t[1]) + Zeta(E, Ee, z - Ee.t[2]))/2));
}

/*
  Parallelograms are represented as triples  [z_0, z_1, z_2],
  corresponding to
  { z_0 + s_1 z_1 + s_2 z_2 | -1/2 <= s_i <= 1/2 }.
*/

/*
  Maximum distance of a point from the centre
  in the parallelogram  R.
*/
par_dist(R) = {
  return(max(abs(R[2] - R[3]), abs(R[2] + R[3]))/2);
}

/*
  Cut the parallelogram  R  into two smaller parallelograms
  along the line through the centre parallel to a shortest side.
*/
par_cut(R) = {
  return(if(abs(R[2]) > abs(R[3]),
            [[R[1] - R[2]/2, R[2]/2, R[3]],
             [R[1] + R[2]/2, R[2]/2, R[3]]],
            [[R[1] - R[3]/2, R[2], R[3]/2],
             [R[1] + R[3]/2, R[2], R[3]/2]]));
}

/*
  Compute the maximum (resp. minimum) of phi on the
  parallelogram  R  to within precision prec.

  The functions below do not check whether the bound  err  is
  valid for the given data.  A correct bound for all cases is

  err = 2 * Ee.W_bound * par_dist(R)
*/
phi_max_recursion(E, Ee, R, prec, z_max, phi_max) = {
  local(phi_z = phi(E, Ee, R[1]),
        W_bound = abs(W(E, Ee, R[1])) + Ee.M * par_dist(R),
        err = 2 * W_bound * par_dist(R));
  if(phi_z > phi_max,
     phi_max = phi_z;
     z_max = R[1]);
  phi_max = max(phi_max, phi_z);
  if(phi_z + err < phi_max + prec,
     return([z_max, phi_max]));
  rec = apply(S -> phi_max_recursion(E, Ee, S, prec, z_max, phi_max),
              par_cut(R));
  return(if(rec[1][2] > rec[2][2], rec[1], rec[2]));
}

phi_min_recursion(E, Ee, R, prec, z_min, phi_min) = {
  local(phi_z = phi(E, Ee, R[1]),
        W_bound = abs(W(E, Ee, R[1])) + Ee.M * par_dist(R),
        err = 2 * W_bound * par_dist(R));
  if(phi_z < phi_min,
     phi_min = phi_z;
     z_min = R[1]);
  if(phi_z - err > phi_min - prec,
     return([z_min, phi_min]));
  rec = apply(S -> phi_min_recursion(E, Ee, S, prec, z_min, phi_min),
              par_cut(R));
  return(if(rec[1][2] < rec[2][2], rec[1], rec[2]));
}

phi_maximum(E, prec) = {
  local(Ee = ellinit_extra(E),
        R0 = [0, E.omega[1], E.omega[2]],
        phi0 = phi(E, Ee, 0));
  return(phi_max_recursion(E, Ee, R0, prec, 0, phi0));
}

phi_minimum(E, prec) = {
  local(Ee = ellinit_extra(E),
        R0 = [0, E.omega[1], E.omega[2]],
        phi0 = phi(E, Ee, 0));
  return(phi_min_recursion(E, Ee, R0, prec, 0, phi0));
}

phi_extrema(E, prec) = {
  return([phi_minimum(E, prec),
          phi_maximum(E, prec)]);
}