\\ The height-pairing work goes wrong if the precision is too small
\p 60

\\ Important Utilities

isratsqr(U)=issquare(numerator(U))&&issquare(denominator(U));
ratsqrt(U)=sqrtint(numerator(U))/sqrtint(denominator(U));
ellas(E)=[E[1],E[2],E[3],E[4],E[5]];
vecge(P,Q)=for(i=1,length(P),if(P[i]<Q[i],return(0)));return(1);

ptslst(P)=
{
\\ List rational points on y^2=P(x)
  local(V,i,j,k);
  V=[];
  for(i=1,25,for(j=-100,100,
    if(gcd(i,j)==1,k=subst(P,x,j/i);if(isratsqr(k),V=concat(V,[[j/i,ratsqrt(k)]])))));
  V;
}

{
ellliset(E,P)=
  local(i,T,P2);P2=[];
  for(i=1,length(P),if(ellheight(E,P[i])>0,P2=concat(P2,[P[i]])));
  T=[P2[1]];
  for(i=2,length(P2),
    if(matdet(ellheightmatrix(E,concat(T,[P2[i]])))>0.001,
      T=concat(T,[P2[i]])));
  T;
}

\\ mestre auxiliary functions

{
mestre_octsqrtres(P)=
\\ compute Q, R such that P=Q^2+R
  local(C,Q,R);
  C=vector(4);
  C[4] = polcoeff(P,7)/2;
  C[3] = (polcoeff(P,6)-C[4]^2)/2;
  C[2] = (polcoeff(P,5)-2*C[3]*C[4])/2;
  C[1] = (polcoeff(P,4)-2*C[2]*C[4]-C[3]^2)/2;
  Q=C[1]+C[2]*x+C[3]*x^2+C[4]*x^3+x^4;
  R=Q^2-P;
}

{
mestre_minimisecubic(R,pts)=
\\ convert y^2=R(x) to a minimal model, transferring the points
  local(lc);
  lc=polcoeff(R,3);
  if(lc!=0,
\\ This is the transformation to a monic cubic
    E=ellinit([0,polcoeff(R,2),0,polcoeff(R,1)*polcoeff(R,3),polcoeff(R,0)*polcoeff(R,3)^2]);
    for(i=1,length(pts),
      pts[i][1]=pts[i][1]*lc;
      pts[i][2]=pts[i][2]*lc);
\\ Now reduce it ...
    eg=ellglobalred(E);
    E=ellchangecurve(E,eg[2]);
    F=ellchangepoint(pts,eg[2]);
    return([E,F]);
  )
}

\\ mestre main functions

{
mestreA(V)=
  local(P,C,Q,R,E,F,eg);
  P=prod(i=1,8,x-V[i]);
  R=mestre_octsqrtres(P);
  print(R);
  T=vector(8,i,[V[i],ratsqrt(subst(R,x,V[i]))]);
  print(T);
  if(polcoeff(R,3)!=0,
    Q=mestre_minimisecubic(R,T);
    F=ellliset(Q[1],Q[2]);
    [length(F),eg[1],ellas(Q[1]),F],[0,0,0,0]);
}

{
mestreB(V)=
  local(P,C,Q,R,E,F,eg);
  P=prod(i=1,8,x-V[i]);
  R=mestre_octsqrtres(P);
  print(R);
  if(polcoeff(R,3)!=0,
    T=ptslst(R);
    Q=mestre_minimisecubic(R,T);
	F=ellliset(Q[1],Q[2]);
	[length(F),ellglobalred(Q[1])[1],ellas(Q[1]),F],[0,0,0,0]);
}

\\ and a running routine

{
mestre_exhaust(bnd)=
  local(Z,ZZ,V,L,M,N);

  Z=vector(25,i,10^100);
  ZZ=vector(25);
  V=vector(6,i,[1,bnd-1]);

  forvec(L=V,
    M=concat(0,L);M=concat(M,bnd);print(M);
    N=mestreB(M);
    if(N[1]>0,
      if(N[2]<Z[N[1]],
      print(M,"\t",N);
      Z[N[1]]=N[2];
      ZZ[N[1]]=N[3];
    ));
  ,2);
  [Z,ZZ];
}