# Maple code to reduce y^2=(quartic in x),  with a rational point,
# to a cubic in two variables. There is an accompanying script to
# further reduce this cubic to canonical form if needed.
# Taken straight from ch. 8 of Cassel's book.
# rusin@math.niu.edu  3/28/96
#
#** NEW  -- 1998/07/20 --
#** An alternative method, not passing to a general cubic, is appended at
#** the end; this seems to work more robustly and efficiently.
#
# Enter quartic...
Q:=1-x^4;
# ... and a rational point on  y^2=Q(x) (i.e.  b0^2=Q(a0).)
a0:=1;
b0:=0;
# We use that rational point only to ensure that  Q  starts with a square:
# indeed if  xbar=x-a0, then Q(x) = Q(xbar + a0) has  b0^2  as its constant
# term when expressed as a polynomial in  xbar, so  b0^2  will be the lead
# coefficient in what follows: if
# ya:=y/(b0*(x-a0)^2);
# xa:=1/(x-a0);
# then  ya^2 = Q(a0+1/xa) * xa^4/Q(a0) = xa^4 + ...
Q1:=simplify(subs(x=a0+1/xa,Q)*xa^4 / b0^2);
# Now attempt to complete the square in  Q1:
G:=xa^2 + (1/2)*coeff(Q1,xa,3)*xa + (1/2)*(coeff(Q1,xa,2) - (1/4)*coeff(Q1,xa,3)^2);
H:=simplify(Q1-G^2);
# So the curve is  ya^2 = G(xa)^2 + H(xa). Factor the difference of squares:
# let T:=ya+G;
# Then ya-G(xa) = H(xa)/T; subtracting shows  2G(xa) = T - H(xa)/T, or
# T^3 - H(xa)*T = 2G(xa)*T^2. We pull out powers of  xa, writing
# S:=xa*T;
fnew:=simplify(T^2*subs(xa=S/T, T-H/T-2*G));
# This is the cubic to solve for  S  and T.
# These are the substitutions made
newT:=simplify(y/(b0*(x-a0)^2)+subs(xa=1/(x-a0),G));
newS:=newT/(x-a0);
# which we can invert if we wish
xorig:=rhs(solve({newT=T, newS=S},{x,y})[1]);
yorig:=rhs(solve({newT=T, newS=S},{x,y})[2]);
#we double check:
factor(simplify(subs({x=xorig,y=yorig},y^2-Q)));
#should give the cubic  fnew  as a factor; conversely
factor(simplify(subs({S=newS,T=newT},fnew)));
# should have  y^2-Q  as a factor.
#
# Now calling in cubic canonical-form reducer
#
f:=subs({S=x,T=y},fnew);
x0:=0;
y0:=0;
read `elliptic.maple`;
#Note that you will need to adjust and rerun that file -- read it.
xall:=simplify(subs({S=xold,T=yold},xorig));
yall:=simplify(subs({S=xold,T=yold},yorig));
#Thus  xall  and  yall  express the original variables  x,y  in the
#quartic  y^2 = Q(x)  as rational functions of the variables  xxx,yyy
#of the canonical cubic, which is  0=
f10;
#now consider running mwrank to find some points on the curve  f10.
#get points on original quartic to get e.g.
subs({xxx=2,yyy=3},{xall,yall});
subs({xxx=-2,yyy=1},{xall,yall});
#
#
==============================================================================

Here's an alternative approach which is perhaps more direct.

Suppose (a0, b0) is a rational point on  y^2=Q(x). Translate to the
origin to suppose  a0=0.

If  b0=0 too  (that is, Q has a rational root)  simply let y=V/U^2
and x=1/U; clearing denominators gives  V^2 = Cubic in U already.

Otherwise, use a transformation of the type
	x = U/V	 y = (a V^2+b UV + c U^2 + d U^3)/V^2
which is invertible if  d <> 0: the inverse is  U = V x  where
V = (Y-a-bx-cx^2)/(dx^3). We need only choose appropriate coefficients
a,b,c,d.

Substitute and clear denominators (V^4); a polynomial of degree 6 results.
Choose  a  to make this polynomial a multiple of  U  (requires
a = +- b0 ); choose b and c to make it a multiple of  U^2  and  U^3  
respectively. Divide by  U^3; what remains is 
	(quadratic in  U  and  V ) = cubic in U.
Choose  d  to make lead coefficients match.
Complete the square (in V) if desired.

Here is a MAPLE routine to do this.


#Given a quartic  QUARTIC in  X  such that  X=a0  is a perfect square,
#return [substitution for X, substitution for Y, resulting elliptic curve] .
#We require global variables  'U' and 'V' to be free so we can report
#the answer!

fixquart:=proc(QUARTIC,X,a0)
global U,V:
local X2,Y,a,b,c,d,W,V0,quart,soln,neweq,allsubs:
#
if U <> 'U' or V <> 'V' then print(`Variables U or V in use! Aborting`):
[]: #this will be the output.
#
elif simplify(subs(X=a0,QUARTIC))=0 then 
print(`Given X coordinate is a root of QUARTIC (Y=0). Thanks!`):
[a0+1/U,V/U^2, 
	V^2=convert(factor(taylor(U^4*subs(X=a0+1/U,QUARTIC),U,5)),polynom)]:
#That's the case  QUARTIC  has a rational _root_ at a0.
#
else
#
quart:=subs(X=X2+a0,QUARTIC-Y^2):
#
subs({X2=U/W,Y=(a*W^2+b*U*W+c*U^2+d*U^3)/W^2},quart*W^4):
#
solve(coeff(",U,0),a); soln:={a="[1]}; factor(subs(soln,""")/U):
#				There should have been two solutions, +- b0
solve(coeff(",U,0),b); soln:={op(soln),b="}; factor(subs(soln,""")/U):
solve(coeff(",U,0),c); soln:={op(soln),c="}; factor(subs(soln,""")/U):
solve(factor((coeff(",U,3)+coeff(",W,2))/d),d);soln:={op(soln),d="};
	factor(subs(soln,""")):  #another solution is  d=0 !
expand(numer(")):"/coeff(",W,2):
#Now you've got a cubic which is only quadratic in  W and lead coeff is 1.
V0:=coeff(expand("),W,1)/2:
subs(W=V-V0,""):#Now it should be V^2 - monic cubic in  U
neweq:=V^2=convert(factor(taylor(V^2-",U,4)),polynom);
#Now you've got  V^2 = cubic in  U
#
subs(soln,{X2=U/W,Y=(a*W^2+b*U*W+c*U^2+d*U^3)/W^2}):
subs(W=V-V0,"):
{X=a0+subs(",X2),Y=subs(",Y)}:
allsubs:=factor("):
#
#check:
print(`Checking...`):
factor(subs(allsubs,QUARTIC-Y^2)):simplify(",{neweq});  
if " <> 0 then print(`Uh-oh -- something went wrong!`): fi:
#
[subs(allsubs,X),subs(allsubs,Y),neweq]: #The answer to be passed.
#
fi: #We checked for  U and V  to be free, remember?
end:

#==============================================================================
#You might want to prettify this answer (There is a general procedure, not
#coded here, which looks for a minimal form for a Weierstrass presentation
#of an elliptic curve. See  APECS.

neweq:="[3]:
#
#We can make it integral:
factor(lcm(denom(coeff(rhs(neweq),U,2\                                      
  ))^6,denom(coeff(rhs(neweq),U,1))^3,denom(coeff(rhs(neweq),U,0))^2));
#Use  ifactor  if these coefficients are integers , not polynomials
#let  R1  be the smallest possible with this last dividing  R1^12;
#(Sorry, I don't know how to do that automagically in Maplese!
lprint(");
R1:=...
nneweq:= V3^2=  # here  V=V3/R1^3
  convert(factor(taylor(subs({U=U3/R1^2},neweq*R1^6),U3,4)),polynom);
#
#Likewise, we can reduce common factors away:
factor(gcd(numer(coeff(rhs(nneweq),U3,2))^6,numer(coeff(rhs(nneweq),U3,1))^3)):
	factor(gcd(",numer(coeff(rhs(nneweq),U3,0))^2));
lprint(");
R2:=... #as before: now the _largest_ R with  R^12 dividing this.
nnneweq:=factor(subs({V3=V4*R2^3,U3=U4*R2^2},nneweq/R2^6));
n4eweq:=V4^2=convert(factor(taylor(rhs(nnneweq),U4,4)),polynom);
#==============================================================================


# Here's the answer as a straightforward formula, for quartics of the
# special form  -y^2  +  A^2*x^4+a3*x^3+a2*x^2+a1*x+a0   for nonzero  a3.

su := 
{x = -(X*Y+a1*X*A+2*a3*a0*A)/(-a3*Y+2*a2*X*A+a1*a3*A+2*X^2*A), y = -(-a3*Y*X^3
-4*a1*X*A^3*a3*a0-4*a0^2*a3^2*A^3-a1^2*X^2*A^3-Y^2*X^2*A-a3^3*Y*a0+2*a2^2*X^3*
A+4*a2*X^4*A+a1^2*a3^2*X*A+3*a1*a3*X^3*A+a1*a3^3*a0*A+2*X^2*a0*a3^2*A-a3*Y*a2*
X^2-a3^2*Y*a1*X+3*a2*X^2*a1*a3*A+2*a2*X*a0*a3^2*A-2*X^2*A^2*Y*a1-4*X*A^2*a3*Y*
a0+2*X^5*A)/(-a3*Y+2*a2*X*A+a1*a3*A+2*X^2*A)^2}:

su_inv:= { X = x*a3-2*y*A+2*A^2*x^2,
Y = -(x^2*a2*a3-2*x*a2*y*A+2*x^3*a2*A^2+x*a1*a3+x^3*a3^2-4*x^2*a3*y*A+4*x^4*
a3*A^2+4*x*y^2*A^2-8*y*A^3*x^3+4*A^4*x^5-A*y*a1+a1*A^2*x^2+a3*a0)/(-y+A*x^2) }:

# The relevant factor from substituting  su  into the quartic is simply

(a0*a3^2-4*a0*a2*A^2+a1^2*A^2)+(-4*a0*A^2+a1*a3)*X+a2*X^2+X^3 - Y^2 ;

# We could of course substitute X = X1 - a2/3  to get an equation without
# quadratic term; the constant term would be
b0:=a0*a3^2-8/3*a0*a2*A^2+a1^2*A^2+2/27*a2^3-1/3*a1*a3*a2 :
# and the coefficient of  X1  would be
b1:= a1*a3-1/3*a2^2-4*a0*A^2 :