Re: a hopeful request

[Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>]

On Sat, Mar 02, 2024 at 09:19:14AM +0000, John Cremona wrote:
> Randall, the formulas for this are in my book. I'm sure Bill must have
> implemented them in libpari, but they are possibly not yet exposed as GP
> functions.

Indeed, I have a GP script implementing the formula from your book.
There is only one thing that I am not sure how to do best:
the reverse (one-to-many) mapping from the curve to the quartic.

Cheers,
Bill

red(q)=q/core(content(q),1)[2]^2
quartic_IJ(g)=
{
  my([a,b,c,d,e]=Vec(g));
  my(ae = a*e, bd = b*d, c2 = c^2, d2=d^2, b2 = b^2);
  my(iI = 12*ae - 3*bd + c2);
  my(iJ = (72*ae  + 9*bd - 2*c2)*c  - 27*a*d2 - 27*b2*e);
  [iI,iJ];
}
IJ_to_ell(V)= ellinit(-27*V);

quartic_to_ell(g)=IJ_to_ell(quartic_IJ(g));

quartic_disc(g)=
{
  my([iI,iJ]=quartic_IJ(g));
  4*iI^3-iJ^2
}
quartic_hessian(g)=
{
  my([A,B,C,D,E]=Vec(g));
  Pol([(3*B^2 - 8*A*C), (4*B*C - 24*A*D), (4*C^2 - 6*B*D - 48*A*E), (4*C*D - 24*B*E), 3*D^2 - 8*C*E]);
}
quartic_GH(g)=
{
  my([iI,iJ]=quartic_IJ(g));
  my(hdd = quartic_hessian(g)'');
  my(Gdd = (4*phi*g'' + hdd)/3);
  my(H = Gdd/12+2/9*(iI-phi^2));
  [G,H]
}

quartic_bracket(g1,g2)=
{
  my(h1=quartic_hessian(g1),h2=quartic_hessian(g2));
  g1*subst(h2,x,y)-subst(g2,x,y)*h1;
}

quartic_isom(g1, g2)=
{
  my(b = quartic_bracket(g1,g2));
  my(S=select(f->poldegree(f)==1,factor(b)[,1]));
  if (#S==0, return(0));
  my([a,b,c,d]=concat([Vecrev(simplify(x),2)|x<-Vecrev(S[1],2)]));
  [-c,-a;d,b];
}

app(Q,m)=
{
  my([u,M]=m);
  subst(Q,x,(M[1,1]*x+M[1,2])/(M[2,1]*x+M[2,2]))*(M[2,1]*x+M[2,2])^4/u^2
}

quartic_pol(g)=
{
  my([iI,iJ]=quartic_IJ(g));
  x^3-3*iI*x+iJ;
}

quartic_cubic(g,phi)=
{
  my([iI,iJ]=quartic_IJ(g));
  my([a,b,c,d,e]=Vec(g));
  (4*a*phi+3*b^2-8*a*c)/3*Mod(1,phi^3-3*iI*phi+iJ);
}

elltocover(C,P)=
{
  my([Q,F]=C,S,T,A);
  F = subst([F[1]*'y^2,F[2]*'y^3],'x,'X); Q=subst(Q,'x,'X);
  S = simplify((F[1]-P[1]*Q)%('y^2-Q));
  T = simplify(F[2]%('y^2-Q));
  A=apply(x->[x,if(1,my(R=subst(T-P[2]*'y*Q,'X,x));-polcoeff(R,0)/polcoeff(R,1))],nfroots(,S));
  if (poldegree(S)<4, my(R=T-P[2]*'y*Q, Y = -pollead(polcoeff(R,0))/pollead(polcoeff(R,1)));
    A=concat(A,[[1,Y,0]]));
  A;
}