/*** PROGRAM velu.gp ***/

\\ Given an elliptic curve E and a torsion point P on E, this
\\ program returns the elliptic curve E' = E/<P> and the isogeny
\\ phi : E --> E' given by Vélu's formulas.
\\
\\ Reference : J. Vélu, Isogénies entre courbes elliptiques,
\\ C. R. Acad. Sc. Paris, Série A, t. 273 (1971), pp. 238-241.

              /*** Example of use  ***/

\\ Let us compute the quotient of E = 11A3 : y^2 + y = x^3 - x^2
\\ by the point P=(0,0), which has order 5 :
\\
\\ e=ellinit([0,-1,1,0,0]);P=[0,0];velu(e,P)
\\
\\ WARNING : make sure that the variables u and v are not used before running the
\\ program. The program first kills u and v, then uses them to describe the isogeny.
\\
\\ The command velu(e,P) returns a vector [eq,phi] where :
\\ 1) eq is a vector [a'1,a'2,a'3,a'4,a'6] giving a Weierstrass equation of E'.
\\ 2) phi is a vector [F(u),G(u,v)] of rational functions, where u and v are
\\     the affine coordinates on E. This means that phi(u,v)=(F(u),G(u,v)).
\\
\\ One can also supply an optional argument d in the function, by typing velu(e,P,d).
\\ The variable d is meant to be the order of P. If d is not specified, the program sets
\\ d=ellorder(e,P). It is useful to specify d when E or P is defined over a number field,
\\ because the function ellorder currently only works for rational points.
\\
\\ The function ellquotient(e,P) initializes the quotient elliptic curve (it is an alias
\\ for ellinit(velu(e,P)[1])). If the function ellidentify is implemented in your version
\\ of PARI, you can get the name of the curve in Cremona's tables and a global minimal
\\ Weierstrass equation by typing ellidentify(ellquotient(e,P)).
\\
\\ Remark : Vélu's isogeny always satisfies phi^* \omega_E' = \omega_E, where
\\ \omega_E and \omega_E' are the differential forms associated to the Weierstrass
\\ equations of E and E'.
\\
\\ To be improved :
\\ 1) When E and/or P is defined over a number field, the program doesn't seem to always work
\\    (the problem lies in Mod(U,eq)). Use velu_coeffs if you only want the equation of the
\\    quotient curve.
\\ 2) Consider non cyclic subgroups.

{
velu(e,P,deg=0)=
local(d,eq,C,U,V,t,w,c);
if(!ellisoncurve(e,P),error("The point is not on the elliptic curve"));
kill(u);
kill(v);
v;u;
d=if(deg,deg,ellorder(e,P));
eq=v^2+e.a1*u*v+e.a3*v-u^3-e.a2*u^2-e.a4*u-e.a6;
C=matrix(d-1,2);
C[1,]=P;
for(k=2,d-1,C[k,]=elladd(e,C[k-1,],P));
U=u+sum(k=1,d-1,ellsub(e,[u,v],C[k,])[1]-C[k,1]);
U=lift(Mod(U,eq));
V=v+sum(k=1,d-1,ellsub(e,[u,v],C[k,])[2]-C[k,2]);
V=lift(Mod(V,eq));
t=sum(k=1,(d-1)/2,6*C[k,1]^2+e.b2*C[k,1]+e.b4);
if(d%2==0,t+=3*C[d/2,1]^2+2*e.a2*C[d/2,1]+e.a4-e.a1*C[d/2,2]);
w=sum(k=1,(d-1)/2,10*C[k,1]^3+2*e.b2*C[k,1]^2+3*e.b4*C[k,1]+e.b6);
if(d%2==0,w+=7*C[d/2,1]^3+(e.b2+2*e.a2)*C[d/2,1]^2+(2*e.b4+e.a4)*C[d/2,1]+e.b6-e.a1*C[d/2,1]*C[d/2,2]);
c=deriv(U,u)*(2*v+e.a1*u+e.a3)/(2*V+e.a1*U+e.a3);
[[e.a1,e.a2,e.a3,e.a4-5*t,e.a6-e.b2*t-7*w],[U,V]]
}

{
ellquotient(e,P)=
ellinit(velu(e,P)[1])
}

{
velu_coeffs(e,P,deg=0)=
local(d,C,t,w,c);
if(!ellisoncurve(e,P),error("The point is not on the elliptic curve"));
d=if(deg,deg,ellorder(e,P));
C=matrix(d-1,2);
C[1,]=P;
for(k=2,d-1,C[k,]=elladd(e,C[k-1,],P));
t=sum(k=1,(d-1)/2,6*C[k,1]^2+e.b2*C[k,1]+e.b4);
if(d%2==0,t+=3*C[d/2,1]^2+2*e.a2*C[d/2,1]+e.a4-e.a1*C[d/2,2]);
w=sum(k=1,(d-1)/2,10*C[k,1]^3+2*e.b2*C[k,1]^2+3*e.b4*C[k,1]+e.b6);
if(d%2==0,w+=7*C[d/2,1]^3+(e.b2+2*e.a2)*C[d/2,1]^2+(2*e.b4+e.a4)*C[d/2,1]+e.b6-e.a1*C[d/2,1]*C[d/2,2]);
[e.a1,e.a2,e.a3,e.a4-5*t,e.a6-e.b2*t-7*w]
}