# Maple routine to put a typical cubic in 2 variables into normal form.
# Assume there is a rational point -- see Cassel's book for picture etc.
# Enter the curve and the point below and tell unix    maple < elliptic.maple
# Or tell maple what  f, x0, y0  are and then          read `elliptic.maple`;
#
# This runs automatically except at the end: run once till it pauses; read
# what vv is, take uu to be its 12-th root (sort of), enter here, and rerun.
# Also note that a few things will make it require human intervention,
# e.g., if the tangent line is already vertical here at the beginning.
#
# rusin@math.niu.edu,  3/9/96
# 
# Begin with a cubic  f(x,y)=0...
#UNCOMMENTME f:=x^3+5*x^2*y+6*y^2*x+y^3-1;
# ... and give the coordinates of a point on it
#UNCOMMENTME x0:=0;
#UNCOMMENTME y0:=1;
# Find the tangent line here
f1:=simplify(subs({x=x0+t,y=y0+m*t},f));
# choose  m  to make t=0 a double root
m1:=solve(coeff(f1,t,1)=0,m);
# now rewrite the tangent line...
f2:=simplify(subs(m=m1,f1));
# ...and find the other point of intersection
t1:=solve(simplify(f2/t^2)=0,t);
x1:=x0+t1;
y1:=y0+t1*m1;
# now cast this new point as the origin
x:=x1+u;
y:=y1+v;
f3:=simplify(f);
# so the original point has u=-t1,v=-t1*m1. Rotate coordinates:
u:=m1*a+b;
v:=-a+m1*b;
f4:=simplify(f3);
# so that the original point has  a=0,b=-t1. Now this line is vertical so
# the equation has to read  b(b+t1)^2+a*Quadratic(a,b)=0.
# Set b=a*c and divide by  a: c*(a*c+t1)^2 + Quadratic(a,c*a), quadratic in a
b:=c*a;
f5:=simplify(f4/a);
# let's complete the square
co2:=coeff(f5,a,2);
co1:=coeff(f5,a,1);
co0:=coeff(f5,a,0);
# so (2 co2 a)^2 + 4 co1 co2 a + 4 co2 co0 = 0, i.e. (2 co2 a + co1)^2=...
a:=(newy-co1)/(2*co2);
f6:=simplify(f5*4*co2);
#Now simplify the cubic
r3:=coeff(f6,c,3);
r2:=coeff(f6,c,2);
c:=d-r2/(3*r3);
f7:=simplify(f6);
# Clear lead coefficient to -1
newy:=-lasty/r3;
d:=-lastx/r3;
f8:=simplify(f7*r3^2);
# make other coefficients integral
s1:=coeff(-f8,lastx,1);
s2:=coeff(lasty^2-f8,lastx,0);
s3:=lcm(denom(s1),denom(s2));
lasty:=yy/s3^3;
lastx:=xx/s3^2;
f9:=simplify(f8*s3^6);
# reduce coefficients by common factors
ss1:=coeff(-f9,xx,1):
ss0:=coeff(yy^2-f9,xx,0):
vv:=ifactor(gcd(ss1^3,ss0^2));
#
#pause(see code for instructions));
#
# When program pauses, quit and rewrite the next uncommented line setting
# uu to the largest integer with  uu^12 dividing  vv. 
# (Thus  uu^4  divides  ss1  and  uu^6 divides  ss0)
# Then comment out the 'pause' line and rerun.
#UNCOMMENTME   uu:=1;
uu:=2^2*3^20*7^10*13*661;
;
yy:=yyy*uu^3:
xx:=xxx*uu^2:
f10:=simplify(f9/uu^6);
#relate to original variables:
xold:=simplify(x);
yold:=simplify(y);