# Last modified 6/14/2008
ESC:=proc () print(`Elliptic Surface Calculator v.1.01. c.1990-2008`) 
end:

interface(typesetting=standard):

`«parameter»`:='t':  # default variable

# These macros are for the readabilitiy of the codes.
macro(
a=`«a»`,
b=`«b»`,
c=`«c»`,
j=`«j»`,#jay=j,
Disc=`«Disc»`,
pg=`«pg»`,
t=`«parameter»`,
curve=`«curve»`
):

# Global variables
# Bad_fibers
# Curve_data

readlib(factors):
readlib(ifactors):
readlib(forget):

Set_parameter:=proc (param_name)
	global t;
	description `Assign var_name as the principal variable.`;
	t:= param_name;
end:

Set_variable:=`Set_parameter`:

elliptic_surface:=proc(a1,a2,a3,a4,a6) 
	global a, b, c, Disc, j, pg, Bad_fibers;
	a[1]:=a1; a[2]:=a2; a[3]:=a3; a[4]:=a4; a[6]:=a6;
        b[2]:=factor(a[1]^2 + 4*a[2]): 
		b[4]:=factor(a[1]*a[3] + 2*a[4]):
        b[6]:=factor(a[3]^2 + 4*a[6]):
        b[8]:=factor(a[1]^2*a[6] + 4*a[2]*a[6] - a[1]*a[3]*a[4] + a[2]*a[3]^2 - a[4]^2):
        c[4]:=factor(b[2]^2 - 24*b[4]):
        c[6]:=factor(-b[2]^3 + 36*b[2]*b[4] - 216*b[6]):
 	Disc:=factor((c[4]^3 - c[6]^2)/1728):
	j:=factor(c[4]^3/Disc):
	if j = 0 then
		pg:=ceil(degree(Disc,t)/12) - 1
	else
		pg:=max(ceil(degree(c[4],t)/4),ceil(degree(Disc,t)/12)) - 1
	fi:
	Bad_fibers:='not_yet':
 	forget(`&++`):forget(`&--`):forget(Add_pts):
 	forget(`&plus`):forget(`&minus`):forget(`&neg`):forget(`&times`):
 	forget(Local_height):forget(Height):forget(Height_pairing):
 print(`\n     The curve you have entered is:`);
 print(y^2 + a[1] * x * y + a[3] * y = sort(x^3 + a[2] * x^2 + a[4] * x + a[6],x));
 print();
end:

Elliptic_surface:=proc()
	if nargs < 5 then
		if degree(args[1],{X,Y}) > 0 then
			elliptic_surface(Extract_coefficients(args[1]));
		elif degree(args[1],{x,y}) > 0 then
			elliptic_surface(Extract_coefficients(subs({x=X,y=Y,z=Z},args[1])));
		fi:
	elif nargs = 5 then
		elliptic_surface(args);
	fi:
end:

Elliptic_curve:=`Elliptic_surface`:
Ein:=`Elliptic_surface`:


Show_data:=proc()
	global a, b, c, Disc, j, pg, Bad_fibers;
	print(y^2 + a[1] * x * y + a[3] * y 
			= sort(x^3 + a[2] * x^2 + a[4] * x + a[6],x));
	print('Discriminant' = Disc);
	print('jay-invariant' = j);

	if Bad_fibers='not_yet' then
		Calculate_bad_fibers();
	fi;

	if pg = 0 then
		print(`\n`);
		print(sprintf(`This is a rational elliptic surface; Oguiso-Shioda type No.%d.`, 
		Oguiso_Shioda_type()));
	elif pg = 1 then
		print(`\n`);
		print(`This is an elliptic K3 surface.`);
	else NULL
	fi;
	
	print(`\n`);
	Show_bad_fibers():
	print(`\n`);
	
	if pg = 0 then
		print(sprintf(
		`The rank of the Mordell-Weil group over C is %d.`,
				Shioda_limit()));
	else
		print(sprintf(
		`The rank of the Mordell-Weil group over C is no greater than %d.`,
				Shioda_limit()));
	fi;
	print(` `);
end:



################################################################
#    Transfromation from cubic equation to Weierstrass form    #
################################################################

# Assume that the equation is given in terms of (x,y,z), and it has a 
# rational point P=(x_0:y_0:z_0).  

# First, move P to (0:1:0) (This means that there is no y^3 term.)
#
step0:=proc(equation,point) local eq, origin, tr, pt;
	eq:=equation;
	pt:=[x,y,z];
	if nargs <> 2 then 
		origin:=[0,1,0];
		if factor(subs({x=origin[1],y=origin[2],z=origin[3]},eq)) <> 0 then
			ERROR(`[0,1,0] is not on the curve!`);
		fi;
	else 
		origin:=point;
		if factor(subs({x=origin[1],y=origin[2],z=origin[3]},eq)) <> 0 then
			ERROR(`The point given is not on the curve!`);
		fi;
	fi;

	if origin[2]=0 then 
		if origin[1]=0 then 
			eq:=mapfactor(subs({y=z,z=y},eq),{x,y,z});
			tr:={x=x,y=z,z=y};
			RETURN(eq,tr,[x,z,y]);
		else
			eq:=mapfactor(subs({x=y,
								y=x,
								z=z+origin[3]/origin[1]*y},eq),
							{x,y,z});
			tr:=factor(subs({	x=y,
								y=x,
								z=z+origin[3]/origin[1]*y}));
			pt:=factor([pt[2], 
						pt[1], 
						pt[3]-origin[3]/origin[1]*pt[1]]);
			RETURN(eq,tr,pt);
		fi;
	fi;
	eq:=mapfactor(numer(normal(
		subs({x=x+origin[1]/origin[2]*y, z=z+origin[3]/origin[2]*y},eq))),
		{x,y,z});
	tr:=factor(subs({x=x+origin[1]/origin[2]*y,y=y, z=z+origin[3]/origin[2]*y}));
	pt:=factor([pt[1]-origin[1]/origin[2]*pt[2], pt[2], pt[3]-origin[3]/origin[2]*pt[2]]);
	RETURN(eq,tr,pt);
end:

# The next step is to move the tangent line at P=[0,1,0] to z=0.  This
# amounts to eliminate x*y^2 term.

step1:=proc(equation,transformation,point) local eq, var, pt, tr;
	eq:=expand(equation);
	var:=subs(transformation,[x,y,z]);
	pt:=point;
	if coeff(coeff(eq,y,2),x) = 0 then
		RETURN(eq,transformation,point);
	fi;
	if coeff(coeff(eq,y,2),z) = 0 then
		var:=factor(subs({x=z,z=x},var));
		pt:=[pt[3],pt[2],pt[1]];
		eq:=mapfactor(subs({x=z,z=x},eq),{x,y,z});
		tr:={x=var[1],y=var[2],z=var[3]};
		RETURN(eq,tr,pt);
	else
		var:=factor(subs({	z=z-coeff(coeff(eq,y,2),x,1)*x,
							x=coeff(coeff(eq,y,2),z,1)*x},var));
		pt:=factor([pt[1]/coeff(coeff(eq,y,2),z,1),pt[2],
				pt[3]+coeff(coeff(eq,y,2),x,1)*pt[1]/coeff(coeff(eq,y,2),z,1)]);
		eq:=mapfactor(
			subs({	z=z-coeff(coeff(eq,y,2),x,1)*x,
					x=coeff(coeff(eq,y,2),z,1)*x},eq),
			{x,y,z});
		tr:={x=var[1],y=var[2],z=var[3]};
		RETURN(eq,tr,pt);
	fi;
end:

# Next, move the third point of intersection with the line z=0 is
# (1:0:0).  (Eliminate x^3.)  If there's no x^2*y term, the given 
# point is a flex, and the procedure ends here.

step2:=proc(equation,transformation,point) local eq, var, tr, pt;
	eq:=expand(equation);
	var:=subs(transformation,[x,y,z]);
	pt:=point;
	if coeff(eq,x,3) = 0 then
		RETURN(eq,transformation,point);
	fi;
	if coeff(coeff(eq,x,2),y) = 0 then
		RETURN(eq,transformation,point);
	else
		var:=cancel(subs({y=y-coeff(eq,x,3)*x,
							x=coeff(coeff(eq,x,2),y,1)*x},var));
		pt:=cancel([pt[1]/coeff(coeff(eq,x,2),y,1),
			pt[2]+coeff(eq,x,3)*pt[1]/coeff(coeff(eq,x,2),y,1), pt[3]]);
		eq:=mapfactor(
			subs({y=y-coeff(eq,x,3)*x,
					x=coeff(coeff(eq,x,2),y,1)*x},eq),
			{x,y,z});
		tr:={x=var[1],y=var[2],z=var[3]};
		RETURN(eq,tr,pt);
	fi;
end:

# Move the tangent line at (1:0:0) to y=0.  (Eliminate x^2*z term.)

step3:=proc(equation,transformation,point) local eq, var, tr, pt;
	eq:=expand(equation);
	var:=subs(transformation,[x,y,z]);
	pt:=point;
	if  coeff(coeff(eq,x,2),y) = 0 then
		RETURN(eq,transformation,point);
	elif coeff(coeff(eq,x,2),z) = 0 then
		var:=cancel(subs({	y=coeff(coeff(eq,x,2),y,1)*y,	
							z=coeff(coeff(eq,x,2),y,1)*z},var));
		pt:=cancel([pt[1],pt[2]/coeff(coeff(eq,x,2),y,1),
					pt[3]/coeff(coeff(eq,x,2),y,1)]);
		eq:=mapfactor(
			subs({	y=coeff(coeff(eq,x,2),y,1)*y,	
					z=coeff(coeff(eq,x,2),y,1)*z},
				eq/coeff(coeff(eq,x,2),y,1)^2),
	    	{x,y,z});
		tr:={x=var[1],y=var[2],z=var[3]};
		RETURN(eq,tr,pt);
	else
		var:=cancel(subs({y=y-coeff(coeff(eq,x,2),z,1)*z,
							z=coeff(coeff(eq,x,2),y,1)*z},var));
		pt:=cancel([pt[1],
			pt[2]+coeff(coeff(eq,x,2),z,1)*pt[3]/coeff(coeff(eq,x,2),y,1),
			pt[3]/coeff(coeff(eq,x,2),y,1)]);
		eq:=mapfactor(
			subs({y=y-coeff(coeff(eq,x,2),z,1)*z,
					z=coeff(coeff(eq,x,2),y,1)*z},eq)
				/coeff(coeff(eq,x,2),y,1),
			{x,y,z});
		tr:={x=var[1],y=var[2],z=var[3]};
		RETURN(eq,tr,pt);
	fi;
end:

# The next step is the heart of this procedure.  Apply the Cremona
# transformation (x,y,z) -> (yz,yx,z^2).  Its inverse is given by
# (X,Y,Z) -> (YZ,X^2,XZ).

step4:=proc(equation,transformation,point) local eq, var, tr, pt;
	eq:=expand(equation);
	var:=subs(transformation,[x,y,z]);
	pt:=point;
	if coeff(coeff(eq,x,2),y) = 0 then
		var:=cancel(subs({	y=Y,
							x=-X*coeff(coeff(eq,y,2),z,1),
							z=Z*coeff(eq,x,3)*coeff(coeff(eq,y,2),z,1)^2},
							var));
		pt:=cancel([-pt[1]/coeff(coeff(eq,y,2),z,1),pt[2],
					pt[3]/(coeff(eq,x,3)*coeff(coeff(eq,y,2),z,1)^2)]);
		eq:=mapfactor(subs({y=Y,
							x=-X*coeff(coeff(eq,y,2),z,1),
							z=Z*coeff(eq,x,3)*coeff(coeff(eq,y,2),z,1)^2},
							eq/coeff(eq,x,3)/coeff(coeff(eq,y,2),z,1)^3),
						{X,Y,Z});
		tr:={x=var[1],y=var[2],z=var[3]};
		RETURN(eq,tr,pt);
	else
		var:=cancel(subs({	x=-coeff(coeff(eq,y,2),z,1)*Y*Z,
							y=X^2,
							z=-coeff(coeff(eq,y,2),z,1)*X*Z},var));
		pt:=cancel([-coeff(coeff(eq,y,2),z,1)*pt[2]*pt[3],
					-coeff(coeff(eq,y,2),z,1)*pt[1]*pt[2],
					pt[3]^2]);
		eq:=mapfactor(
			subs({	x=-coeff(coeff(eq,y,2),z,1)*Y*Z,
					y=X^2,
					z=-coeff(coeff(eq,y,2),z,1)*X*Z},eq)
				/(coeff(coeff(eq,y,2),z,1)^2*X^2*Z),
			{X,Y,Z});
		tr:={x=var[1],y=var[2],z=var[3]};
		RETURN(eq,tr,pt);
	fi;
end:

# At this point, the equation should be of the form
# Y^2Z + a1 XYZ + a3 YZ^2 - X^3 - a2 X^2Z - a4 XZ^2 - a6 Z^3.

# Eliminate a1 and a3.

step5:=proc(equation,transformation,point) local eq, var, tr, pt, a1, a3;
	eq:=expand(equation);
	a1:=coeff(coeff(coeff(eq,Y,1),X,1),Z,1);
	a3:=coeff(coeff(eq,Y,1),Z,2);
	var:=subs(transformation,[x,y,z]);
	pt:=point;
	if modp(coeff(eq,Y),2) = 0 then 
		var:=cancel(subs({Y=Y-(a1*X+a3*Z)/2},var));
		pt:=cancel([pt[1],pt[2]+(a1*pt[1]+a3*pt[3])/2,pt[3]]);
		eq:=mapfactor(subs({Y=Y-(a1*X+a3*Z)/2},eq),{X,Y,Z});
	else
		var:=cancel(subs({X=2*X,Y=Y-(a1*X+4*a3*Z),Z=8*Z},var));
		pt:=cancel([pt[1]/2,
					pt[2]+(a1*pt[1]+a3*pt[3])/2,
					pt[3]/8]);
		eq:=mapfactor(
			subs({X=2*X,Y=Y-(a1*X+4*a3*Z),Z=8*Z},eq/2^3),
			{X,Y,Z});
	fi;
	tr:={x=var[1],y=var[2],z=var[3]};
	RETURN(eq,tr,pt);
	end:

# Eliminate a2.

step6:=proc(equation,transformation,point) local eq, var, tr, pt, a2;
	eq:=expand(equation);
	a2:=coeff(coeff(eq,X,2),Z,1);
	var:=subs(transformation,[x,y,z]);
	pt:=point;
	if modp(coeff(eq,X^2),3) = 0 then 
		var:=cancel(subs({X=X+a2*Z/3},var));
		pt:=cancel([pt[1]-a2*pt[3]/3,pt[2],pt[3]]);
		eq:=mapfactor(
			subs({X=X+a2*Z/3},eq),
			{X,Y,Z});
	else
		var:=cancel(subs({X=3*X+9*a2*Z,Z=27*Z},var));
		pt:=cancel([pt[1]/3-a2*pt[3]/9,pt[2],pt[3]/27]);
		eq:=mapfactor(
			subs({X=3*X+9*a2*Z,Z=27*Z},eq/3^3),
			{X,Y,Z});
	fi;
	tr:={x=var[1],y=var[2],z=var[3]};
	RETURN(eq,tr,pt);
	end:

# Eliminate extra factors from a4 and a6, i.e.,
# Y^2 Z = X^3 + u^4 a4 X Z^2 + u^6 a6 Z^3
#  ==> Y^2 Z = X^3 + a4 X Z^2 + a6 Z^3.
# For now, you have to input 'u' manually.

readlib(factors):
readlib(ifactors):

step7:=proc(equation,transformation,point) 
	local eq, c4, c6, nc4, nc6, list, u, var, tr, pt;
	c4:=-coeff(coeff(equation,X,1),Z,2);
	c6:=-coeff(equation,Z,3);
	nc4:=op(1,factors(c4))  # numerical part
				/(denom @ `*` @ op @ map)(x->x[1]^x[2],op(2,factors(c4)));
	nc6:=op(1,factors(c6))
				/(denom @ `*` @ op @ map)(x->x[1]^x[2],op(2,factors(c6)));
	list:=map(x->x[1], op(2,factors(gcd(c4,c6)))); # polynomial factors
	if c4 = 0 then
		u:=(`*` @ op @ map)
				((x,c_6)->	numer(x)^floor(ord1(c_6,x)/6),
					list,c6);
		list:=map(x->x[1], 
				op(2,ifactors(nc6)));  # numerical powers
		u:=u * (`*` @ op @ map)
				((x,c_6)->x^floor(ord2(c_6,x)/6),
					list,c6);	
		var:=subs(transformation,[x,y,z]);
		pt:=point;
		var:=cancel(cancel(subs({Y=Y*u^3,X=X*u^2},var)));
		pt:=cancel([pt[1]/u^2,pt[2]/u^3,pt[3]]);
		eq:=Y^2*Z - (X^3 + normal(c6/u^6)*Z^3);
		tr:={x=var[1],y=var[2],z=var[3]};	
		RETURN(eq,tr,pt);
	elif c6 = 0 then  # has not been tested
		u:=(`*` @ op @ map)
				((x,c_4)->	numer(x)^floor(ord1(c_4,x)/4),
					list,c4);
		list:=map(x->x[1], 
				op(2,ifactors(nc4*nc6)));  # numerical powers
		u:=u * (`*` @ op @ map)
				((x,c_4)->x^floor(ord2(c_4,x)/4),
					list,c4,c6);	
		var:=subs(transformation,[x,y,z]);
		pt:=point;
		var:=cancel(cancel(subs({Y=Y*u^3,X=X*u^2},var)));
		pt:=cancel([pt[1]/u^2,pt[2]/u^3,pt[3]]);
		eq:=Y^2*Z - (X^3 + normal(c4/u^4)*X*Z^2);
		tr:={x=var[1],y=var[2],z=var[3]};	
		RETURN(eq,tr,pt);
	else
		u:=(`*` @ op @ map)
				((x,c_4,c_6)->
						numer(x)^min(floor(ord1(c_4,x)/4),floor(ord1(c_6,x)/6)),
					list,c4,c6);
		list:=map(x->x[1], 
				op(2,ifactors(nc4*nc6)));  # numerical powers
		u:=u * (`*` @ op @ map)
				((x,c_4,c_6)->x^min(floor(ord2(c_4,x)/4),floor(ord2(c_6,x)/6)),
					list,c4,c6);	
		var:=subs(transformation,[x,y,z]);
		pt:=point;
		var:=cancel(cancel(subs({Y=Y*u^3,X=X*u^2},var)));
			 # The second "cancel" is to remove numerical coefficients.
		pt:=cancel([pt[1]/u^2,pt[2]/u^3,pt[3]]);
		eq:=Y^2*Z - (X^3 + normal(c4/u^4)*X*Z^2 + normal(c6/u^6)*Z^3);
		tr:={x=var[1],y=var[2],z=var[3]};	
		RETURN(eq,tr,pt);
	fi;
	end:

step7bis:=proc(equation,transformation,point) 
	local eq, c4, c6, list, u, var, tr, pt;
	c4:=-coeff(coeff(equation,X,1),Z,2);
	c6:=-coeff(equation,Z,3);
	list:=map(x->x[1], op(2,factors(gcd(c4,c6)))); # polynomial factors
	u:=(`*` @ op @ map)
			((x,c_4,c_6)->x^min(floor(ord1(c_4,x)/4),floor(ord1(c_6,x)/6)),
				list,c4,c6);
	pt:=point;
	pt:=cancel([pt[1]/u^2,pt[2]/u^3,pt[3]]);
	eq:=Y^2*Z - (X^3 + normal(c4/u^4)*X*Z^2 + normal(c6/u^6)*Z^3);
	RETURN(eq,pt);
	end:
	
ord1:= proc (p1, p2)
	local list;
	if p1 = 0 then RETURN(infinity) fi;
	list:=op(2,factors(p1)):
	(`+` @ op @ map)((x,y) -> if y = x[1] then x[2] else 0 fi,list,p2);
end:

ord2:= proc (poly, q) # numerical coeffs
  	local list;
  	list:=op(2,ifactors(
				op(1,factors(poly))
				/(denom @ `*` @ op @ map)(x->x[1]^x[2],op(2,factors(poly)))
			)):
  	(`+` @ op @ map)((x,y) -> if y = x[1] then x[2] else 0 fi,list,q);
end:

cancel:=proc (point) local pt, pgcd, ppcm;
	pt:=point;
	ppcm:=lcm(lcm(denom(pt[1]),denom(pt[2])),denom(pt[3]));
	pt:=factor([ppcm*pt[1],ppcm*pt[2],ppcm*pt[3]]);
	pgcd:=gcd(gcd(pt[1],pt[2]),pt[3]);
	pt:=factor([pt[1]/pgcd,pt[2]/pgcd,pt[3]/pgcd]);
	end:

	
# mapfactor:=proc(expression, variables) 
#  		map(factor,collect(
#  			expand(expression) 
#  			- subs(map(x->x=0,variables),simplify(expression)),
#  		variables,'distributed'))
#  		+ factor(subs(map(x->x=0,variables),simplify(expression)))
#  	   end:	
	
Cubic_to_Weierstrass2:=proc(curve,point) 
if nargs=1 then
step6(step5(step4(step3(step2(step1(step0(curve)))))))
else step6(step5(step4(step3(step2(step1(step0(curve,point))))))) fi; 
end:

Cubic_to_Weierstrass1:=proc(curve,point) 
if nargs=1 then
step4(step3(step2(step1(step0(curve)))))
else step4(step3(step2(step1(step0(curve,point))))) fi; 
end:

Cubic_to_Weierstrass:=proc(curve,point) 
if nargs=1 then
step7(step6(step5(step4(step3(step2(step1(step0(curve))))))))
else step7(step6(step5(step4(step3(step2(step1(step0(curve,point)))))))) fi; 
end:

CtoW:=`Cubic_to_Weierstrass`:


Extract_coefficients:=proc(equation,var) local eq, lc;
	if nargs = 1 then
		eq:=expand(subs({Z=1},equation));
		lc:=coeff(eq,Y,2);
		if coeff(eq,X,3) + lc <> 0 then
			ERROR(`Equation is not of the form Y^2Z - X^3 -...`);
		fi;
		eq:=mapfactor(eq/lc,{X,Y});
	else
		eq:=expand(subs({var[1]=X,var[2]=Y,var[3]=1},equation));
		lc:=coeff(eq,Y,2);
		if coeff(eq,X,3) + lc <> 0 then
			ERROR(`Equation is not of the form Y^2Z - X^3 -...`);
		fi;
		eq:=mapfactor(eq/lc,{X,Y});
	fi;
	RETURN(	coeff(coeff(eq,Y,1),X,1),
			-coeff(eq,X,2),
			coeff(coeff(eq,Y,1),X,0),
			-coeff(coeff(eq,X,1),Y,0),
			-coeff(coeff(eq,X,0),Y,0) );
end:

Extract_coeffs:=`Extract_coefficients`:

################################################################
#   Transfromation from  equation to Weierstrass form   #
################################################################

# Quartic equation is assumed to be
# V^2 = A*U^4 + B*U^3 + C*U^2 + D*U + E

Quartic_to_Weierstrass:=proc(equation,point) 
	local quartic, eq, p, q, A, B, C, D, E, tr, tr_inv, a1,a2,a3,a4,a6,b2, b4, b6, list, u, var, pt;
	
	if degree(equation, V) > 0 then
		quartic:=-subs(V=0,equation)/coeff(equation,V,2);
	else
		quartic:=equation;
	fi;
	p:=point[1];
	q:=point[2];
	A:=coeff(quartic,U,4);
	B:=coeff(quartic,U,3);
	C:=coeff(quartic,U,2);
	D:=coeff(quartic,U,1);
	E:=coeff(quartic,U,0);
	
	if q = 0 then
		if factor(rem(quartic,U-p,U)) = 0 then
			A:=coeff(quo(quartic,U-p,U),U,3);
			B:=coeff(quo(quartic,U-p,U),U,2);
			C:=coeff(quo(quartic,U-p,U),U,1);
			D:=coeff(quo(quartic,U-p,U),U,0);
			eq:=mapfactor( 
				Y^2*Z-X^3-(2*B*p+3*A*p^2+C)*X^2*Z
				-(B+3*A*p)*(C*p+D+B*p^2+A*p^3)*X*Z^2
				-A*(C*p+D+B*p^2+A*p^3)^2*Z^3,
				{X,Y,Z});
			b2:=(2*B*p+3*A*p^2+C);
			b4:=(B+3*A*p)*(C*p+D+B*p^2+A*p^3);
			b6:=A*(C*p+D+B*p^2+A*p^3)^2;
			list:=map(x->x[1], op(2,factors(gcd(b2,gcd(b4,b6))))); # polynomial factors
			u:=(`*` @ op @ map)
			        ((x,b_2,b_4,b_6)->
			                x^min(floor(ord1(b_2,x)/2),floor(ord1(b_4,x)/4),floor(ord1(b_6,x)/6)),
				 list,b2,b4,b6);
			eq:=Y^2*Z - (X^3 + factor(b2/u^2)*X^2*Z 
				        + factor(b4/u^4)*X*Z^2 + factor(b6/u^6)*Z^3);
			tr:=simplify(
				{X =factor((C*p+D+B*p^2+A*p^3)*(U-p)/u^2),
			 	Y = factor((C*p+D+B*p^2+A*p^3)*V/u^3),
			 	Z = (-U+p)^2
				});
			tr_inv:=simplify(subs( {X=X/Z,Y=Y/Z},
				{U = ((C*p+D+B*p^2+A*p^3)+p*X)/(X*u^2),
				 V = Y*u^3*(C*p+D+B*p^2+A*p^3)/(X*u^2)^2
				}));			
			RETURN(eq,tr,tr_inv);
		else 
			ERROR(`Given point is not on the curve!`); 
		fi;
	else
		if factor(subs({U=p,V=q}, V^2 - quartic)) = 0
		then
			eq:=mapfactor(
				Y^2*Z + (D+2*C*p+3*B*p^2+4*A*p^3)*X*Y*Z +
				2*(B+4*A*p)*q^4*Y*Z^2 
				- ( X^3 
					+ ((C+3*B*p+6*A*p^2)*q^2
					-1/4*(D+2*C*p+3*B*p^2+4*A*p^3)^2)*X^2*Z
					- 4*A*q^6*X*Z^2
					+ (-4*((C+3*B*p+6*A*p^2)*q^2
					-1/4*(D+2*C*p+3*B*p^2+4*A*p^3)^2)*q^6*A)*Z^3
				  ),
				{X,Y,Z});
			a1:=(D+2*C*p+3*B*p^2+4*A*p^3);
			a2:=((C+3*B*p+6*A*p^2)*q^2
					-1/4*(D+2*C*p+3*B*p^2+4*A*p^3)^2);
			a3:=2*(B+4*A*p)*q^4;
			a4:=- 4*A*q^6;
			a6:=(-4*((C+3*B*p+6*A*p^2)*q^2
					-1/4*(D+2*C*p+3*B*p^2+4*A*p^3)^2)*q^6*A);
			list:=map(x->x[1], op(2,factors(gcd(a1,gcd(a2,gcd(a3,gcd(a4,a6))))))); # polynomial factors
			u:=(`*` @ op @ map)
			        ((x,a_1,a_2,a_3,a_4,a_6)->
			                x^min(floor(ord1(a_1,x)/1),floor(ord1(a_2,x)/2),
			                floor(ord1(a_3,x)/3),floor(ord1(a_4,x)/4),floor(ord1(a_6,x)/6)),
				 list,a1,a2,a3,a4,a6);
			eq:=Y^2*Z +factor(a1/u)*X*Y*Z + factor(a3/u^3)*Y*Z^2 
				        - (X^3 + factor(a2/u^2)*X^2*Z 
				        + factor(a4/u^4)*X*Z^2 + factor(a6/u^6)*Z^3);
			tr:=simplify(
				{X = factor(q^2*(2*q*(V+q)+(D+2*C*p+3*B*p^2+4*A*p^3)*(U-p))*(U-p)/u^2),
				 Y = factor(q^3*(4*q^2*(V+q)+2*q*((D+2*C*p+3*B*p^2+4*A*p^3)*(U-p)
						+(C+3*B*p+6*A*p^2)*(U-p)^2)
						-1/2*(D+2*C*p+3*B*p^2+4*A*p^3)^2*(U-p)^2/q)/u^3),
				 Z = (U-p)^3
				});
			tr_inv:=simplify(subs( {X=X/Z,Y=Y/Z},
				{U = (2*q*((X*u^2)/(q^2)+C+3*B*p+6*A*p^2) 
					-1/2*(D+2*C*p+3*B*p^2+4*A*p^3)^2/q)*q^3/(Y*u^3)+p,
				 V = subs(U = (2*q*((X*u^2)/(q^2)+C+3*B*p+6*A*p^2) 
					-1/2*(D+2*C*p+3*B*p^2+4*A*p^3)^2/q)*q^3/(Y*u^3)+p,
					-q+1/2*(U-p)*((U-p)*(X*u^2)/(q^2)-D-2*C*p-3*B*p^2-4*A*p^3)/q)
				}));			
			RETURN(eq,tr,tr_inv);
		else 
			ERROR(`Given point is not on the curve!`); 
	        fi;
	        fi;
end:

QtoW:=`Quartic_to_Weierstrass`:

Jacobian:=proc(quartic) 
	local A, B, C, D, E, i, J, Jac;
	A:=coeff(quartic,U,4);
	B:=coeff(quartic,U,3);
	C:=coeff(quartic,U,2);
	D:=coeff(quartic,U,1);
	E:=coeff(quartic,U,0);
	i:=-27*(12*A*E-3*B*D+C^2);
	J:=-27*(72*A*C*E+9*B*C*D-27*A*D^2-27*B^2*E-2*C^3);
	Jac:=mapfactor(Y^2*Z - ( X^3 + i*X*Z^2 + J*Z^3), {X,Y});
	RETURN(Jac);
end:
	
	

# Quartic_to_Weierstrass:='readlib('Quartic_to_Weierstrass',``.savelibname.` sources:quartic.maple`)':
# Jacobian:='readlib('Jacobian',``.savelibname.` sources:quartic.maple`)':


################################################################
#                     Group operations                         #
################################################################
`&++`:=proc(P1,P2)  # To do: add option `check if it's on`.
	local slope, x1,x2,y1,y2,x3,y3;
	option remember;
	x1:=P1[1];y1:=P1[2];x2:=P2[1];y2:=P2[2];

	if   P1 = `Origin` then RETURN(P2)
	elif P2 = `Origin` then RETURN(P1)
	fi;
	
	if simplify(x1 - x2) = 0 
	then if simplify(y1 - y2) = 0 and ( 2*y1 + a[1]*x1 + a[3] ) <> 0
	     then slope:=( 3*x1^2 + 2*a[2]*x1 + a[4] - a[1]*y1 )
		 				/( 2*y1 + a[1]*x1 + a[3] );
	     else RETURN('Origin');	
	     fi;
	else slope:=( y2 - y1 )/( x2 - x1 );
	fi;
	x3:=factor(slope^2 + a[1]*slope - a[2] - x1 - x2);
    y3:=factor(-slope*(x3-x1) - y1 - a[1]*x3 - a[3]);
	RETURN([x3,y3]);
end:

`&plus`:=`&++`:

`&neg`:=proc(P)
	if P = `Origin` then RETURN(P)
	fi;
	
	[P[1],-(P[2] + a[1]*P[1] + a[3])]
end:

`&negative`:=`&neg`:

`&--`:=proc(P1,P2)
	option remember;

	P1 &++ (&neg P2)
end:

`&mius`:=`&--`:

`&times`:= proc (int, point)
	local i, result;
	result:= `Origin`;
	for i from 1 to int do
		result:= result &++ point
	od;
	RETURN(result);
end:


Is_on:=proc(P)
	local x, y;
	x:=P[1];y:=P[2];
	if simplify(y^2 + a[1]*x*y + a[3]*y 
					- (x^3 + a[2]*x^2 + a[4]*x + a[6]) ) = 0
	then print(`True`);
	else print(`False`);
	fi;
end:

################################################################
#                         abscis                               #
################################################################
Absc:=proc(x)
 simplify(sqrt(factor(x^3+a[2]*x^2+a[4]*x+a[6]),symbolic),symbolic);
end:

################################################################
#                  Calculation of bad fibers                   #
################################################################

Bad_places := proc ()
	global Bad_fibers;
	local list;
	description 
		`Determine the places of bad fibers.  The result is given \
		in the format [[place,degree in discrim],...].`;
		
	list:=map(x -> if degree(x[1],t) > 0 then x fi,
				op(2,factors(Disc)));
	list:=sort(map(x->[sort(op(1,x)),op(2,x)],list),
			(a,b)->if op(2,a) > op(2,b) then true else false fi);
	if  12*(pg + 1) - degree(Disc,t) <> 0 then
		list := [op(list),[infinity,12*(pg + 1) - degree(Disc,t)]]; 
	fi;
	RETURN(list);
end:


Bad_fiber_type := proc (place,ord_disc) 
	global Bad_fibers;
	local ord_c4;
	description
		` \
		Bad_fiber_type[1] = type (1 = I, 2 = II,...) \
		Bad_fiber_type[2] = star (0 = no star, 1 = with star) \
		Bad_fiber_type[3] = order in discrim. (= Euler number) \
		Bad_fiber_type[4] = Kodaira type in printable form \
		Bad_fiber_type[5] = number of components \
		Bad_fiber_type[6] = number of simple components \
		Bad_fiber_type[7] = group structure ([a,b] = Z/aZ x Z/bZ) \
		Note: 	I[b] -> G_m x Z/bZ  \
				II   -> G_a         \
		 		III  -> G_a x Z/2Z  \
				IV   -> G_a x Z/3Z  \
				I[2b]*   -> G_a x (Z/2Z)^2  \
				I[2b+1]* -> G_a x Z/4Z      \
				II*  -> G_a         \
				III* -> G_a x Z/2Z  \
				IV*  -> G_a x Z/3Z  \
		`;
		
	if place = infinity then
		ord_c4:= 4*(pg + 1) - degree(c[4],t);
	else
		ord_c4:=Ord(c[4],place);
	fi;
	
	if 
		ord_disc = 0  # type I[0] (= good fiber)
		then RETURN(1, 0, 0, I[0], 1, 1,'NA')
	elif
		ord_c4 = 0    # type I[b]
		then RETURN(1, 0, ord_disc, I[ord_disc], ord_disc, ord_disc, [0,ord_disc])
	elif 
		ord_disc > 3*ord_c4  # type I[b]*
		then 
		if modp(ord_disc, 2) = 0 then
			RETURN(1, 1, ord_disc, I[ord_disc-6]^`*`, ord_disc-6+5, 4, [2,2])
		else
			RETURN(1, 1, ord_disc, I[ord_disc-6]^`*`, ord_disc-6+5, 4, [4])
		fi;
	elif
		ord_disc = 2
		then RETURN(2, 0, 2, II[``], 1, 1,[1])
	elif
		ord_disc = 3
		then RETURN(3, 0, 3, III[``], 2, 2,[2])
	elif
		ord_disc = 4
		then RETURN(4, 0, 4, IV[``], 3, 3,[3])
	elif
		ord_disc = 6
		then RETURN(1, 1, 6, I[0]^`*`, 5, 4,[2,2])
	elif
		ord_disc = 8
		then RETURN(4, 1, 8, IV[``]^`*`, 7, 3,[3])
	elif
		ord_disc = 9
		then RETURN(3, 1, 9, III[``]^`*`, 8, 2,[2])
	elif
		ord_disc = 10
		then RETURN(2, 1, 10, II[``]^`*`, 9, 1,[1])
	else
		ERROR(`The given equation is not minimal!`);
	fi;
end:

Calculate_bad_fibers:=proc ()
	global Bad_fibers;
	description 
		`Calculate_bad_fibers() computes the places and the types \
		of bad fibers of the elliptic surface. It sets the global \
		variable 'Bad_fibers' whose format is \
		[[place, degree, type(1,2,...), star(0 or 1), deg Disc, \
		Kodaira symbol,# of comps, # of simple comps],...]. \
		For example [t^2+1,2,1,1,7,I[1]^*,6,4] means that the surface \
		has a bad fiber of type I[1]* at t=i and -i, and each I[1]* fiber \
		contains 6 irreducible components and 4 components with multiplicity \
		one.`;

	Bad_fibers:=
		map(pl->if pl[1] = infinity then
					[pl[1],1,Bad_fiber_type(pl[1],pl[2])] 
				else
					[pl[1],degree(pl[1],t),Bad_fiber_type(pl[1],pl[2])]
				fi,
			Bad_places());
end:
	
Show_bad_fibers:=proc ()
	global Bad_fibers;
	local bad_fiber_table, i, n;
	description
		`Print the list of bad fibers in a nice format.`;
	
	if Bad_fibers='not_yet' then
		Calculate_bad_fibers();
	fi;
	n:=nops(Bad_fibers);
	bad_fiber_table:=array(1..n,1..7);
	for i to n	do 
		bad_fiber_table[i,2]:= Bad_fibers[i][6] od:  # Kodaira type
	for i to n do 
		if Bad_fibers[i][1] = infinity then
			bad_fiber_table[i,6]:= (t = infinity)
		elif Bad_fibers[i][2] = 1 then
			bad_fiber_table[i,6]:= (t = solve(Bad_fibers[i][1]=0,t) )
		else		
			bad_fiber_table[i,6]:= ( numer(Bad_fibers[i][1]) = 0 )
		fi
	od:
	for i to n do bad_fiber_table[i,1]:=`     ` od:
	for i to n do bad_fiber_table[i,3]:=`  ` od:
	for i to n do bad_fiber_table[i,4]:=` at ` od:
	for i to n do bad_fiber_table[i,5]:=`  ` od:
	for i to n do bad_fiber_table[i,7]:=`     ` od:
	print(`Singular fibers` = op(bad_fiber_table));
	RETURN(NULL):
end:


# Bad_fiber_list:=proc ()
# 	global Bad_fibers;
# 	local i;
# 	[seq(Bad_fibers[i,2],i=op(op(Bad_fibers))[2,1])];
# end:


Repetition := proc (poly)
	if poly = infinity then 1 
	elif poly = 0 then infinity
	else degree(poly,t)
	fi;
end:

		 	

Possible_torsion:= proc()  # Bug!!!!
	global Bad_fibers;
	(`*` @ op @ map) (
		a -> if irem( (igcd @ op @ map)
						(ilcm @ op,
					 	 map(Group_structure @ Bad_fiber_type @ op,Bad_places())
						),
					 a[1]) = 0
			 then a[1]^floor(a[2]/2) else 1 fi,
		op(2,(ifactors @ `*` @ op @ map)
				(a-> No_of_simple_components(a[1])^a[2],
				 map(a -> [(Bad_fiber_type @ op)(a), Repetition(a[1])],
				 	 Bad_places()))
			)
	)
end:

################################################################
#      Oguiso-Shioda type for rational elliptic surface        #
################################################################


Shioda_limit := proc ()
	global pg ; 
	
	if Bad_fibers='not_yet' then
		Calculate_bad_fibers();
	fi;

	8 + 10*pg - # \sum (No_of_components - 1)*degree 
	    (Add_all @ map)(bad_fiber -> (bad_fiber[7] - 1)*bad_fiber[2],Bad_fibers);
end:


Oguiso_Shioda_type := proc ()
		global a, b, c, Disc, j, pg, Bad_fibers;
		local rank, detT;
	if pg <> 0 then ERROR(`The surface is not a rational surface.`) fi;
	
	rank:=Shioda_limit();
	detT:=(`*` @ op @ map)  # \prod (number of simple comp.)^degree
				(bad_fiber -> bad_fiber[8]^bad_fiber[2],  Bad_fibers);
					 
	if rank > 6 then RETURN(9 - rank)
	elif rank = 6 then 
		if   detT = 3 then RETURN(3) 
		elif detT = 4 then RETURN(4)
		fi;
	elif rank = 5 then 
		if   detT = 4 then RETURN(5) 
		elif detT = 6 then RETURN(6)
		elif detT = 8 then RETURN(7)
		fi;
	elif rank = 4 then 
		if   detT =  5 then RETURN(8) 
		elif detT =  4 then RETURN(9)
		elif detT =  8 then RETURN(10)
		elif detT =  9 then RETURN(11)
		elif detT = 12 then RETURN(12)
		elif detT = 16 then 
			if nops(op(2,factors(x^3-27*c[4]*x-54*c[6]))) = 2
			then RETURN(13)
			else RETURN(14)
			fi;
		fi;
	elif rank = 3 then 
		if   detT =  6 then RETURN(15) 
		elif detT =  4 then RETURN(16)
		elif detT = 10 then RETURN(17)
		elif detT =  8 then RETURN(18)
		elif detT = 12 then RETURN(19)
		elif detT = 18 then RETURN(20)
		elif detT = 16 then 
			if nops(op(2,factors(x^3-27*c[4]*x-54*c[6]))) = 2
			then RETURN(21)
			else RETURN(22)
			fi;
		elif detT = 24 then RETURN(23)
		elif detT = 32 then RETURN(24)
		fi;
	elif rank = 2 then 
		if   detT =  7 then RETURN(25) 
		elif detT =  4 then RETURN(26)
		elif detT =  3 then RETURN(27)
		elif detT = 12 then 
			if nops(op(2,factors(x^3-27*c[4]*x-54*c[6]))) = 2
			then RETURN(28)
			else if member(I[0]^`*`,map(op,Bad_fibers)) then RETURN(29)
				 else RETURN("A bug fix is under way.  Sorry.") #RETURN(32)
				 fi;
			fi;
		elif detT =  8 then RETURN(30)
		elif detT = 15 then RETURN(31)
		elif detT = 20 then RETURN(33)
		elif detT = 16 then 
			if nops(op(2,factors(x^3-27*c[4]*x-54*c[6]))) = 2
			then if member(I[0]^`*`,map(op,Bad_fibers)) then RETURN(34)
				 else RETURN(35)
				 fi;
			else RETURN(36)
			fi;
		elif detT = 24 then RETURN(37)
		elif detT = 32 then RETURN(38)
		elif detT = 27 then RETURN(39)
		elif detT = 36 then RETURN(40)
		elif detT = 48 then RETURN(41)
		elif detT = 64 then RETURN(42)
		fi;
	elif rank = 1 then 
		if   detT =  2 then RETURN(43) 
		elif detT =  4 then RETURN(46)
		elif detT =  6 then RETURN(49)
		elif detT = 8 then 
			if nops(op(2,factors(x^3-27*c[4]*x-54*c[6]))) = 1
			then RETURN(45)
			else if member(I[0]^`*`,map(op,Bad_fibers)) then RETURN(48)
				 else RETURN(44)
				 fi;
			fi;
		elif detT = 14 then RETURN(47)
		elif detT = 12 then RETURN(50)
		elif detT = 18 then RETURN(51)
		elif detT = 16 then 
			if member(I[0]^`*`,map(op,Bad_fibers)) 
			then RETURN(54)
			else RETURN(52)
			fi;
		elif detT = 24 then RETURN(53)
		elif detT = 20 then RETURN(55)
		elif detT = 30 then RETURN(56)
		elif detT = 32 then
			if member(I[0]^`*`,map(op,Bad_fibers)) 
			then RETURN(57)
			else RETURN(58)
			fi;
		elif detT = 48 then RETURN(59)
		elif detT = 64 then RETURN(60)
		elif detT = 54 then RETURN(61)
		fi;
	elif rank = 0 then 
		if   detT =  1 then RETURN(62) 
		elif detT =  9 then 
			if member(I[9],map(op,Bad_fibers)) 
			then RETURN(63)
			else RETURN(69)
			fi;
		elif detT =  4 then 
			if member(I[4]^`*`,map(op,Bad_fibers)) 
			then RETURN(64)
			else RETURN(65)
			fi;
		elif detT = 36 then RETURN(66)
		elif detT = 25 then RETURN(67)
		elif detT = 81 then RETURN(68)
		elif detT = 16 then 
			if member(I[8],map(op,Bad_fibers)) 
			then RETURN(70)
			elif member(I[2]^`*`,map(op,Bad_fibers)) 
			then RETURN(71)
			elif member(I[1]^`*`,map(op,Bad_fibers)) 
			then RETURN(72)
			else RETURN(73)
			fi;
		elif detT = 64 then RETURN(74)
		fi;
	fi;
end:

################################################################
#                      Height functions                        #
################################################################

Local_height:= proc(point, prime)
	global a, b, c, Disc, j, pg, Bad_fibers;
	local A, B, C, D, X, x, y, a_inf, b_inf, c_inf, x_inf, y_inf, i;
	description 
		`Local_height([x,y],p) computes the local height of the \
		point [x,y] for the place p. It is assumed that the \
		Weierstrass equation is minimal at  p.`;
	option remember;
	
	if prime <> infinity then
		x:=point[1]; y:=point[2];
####################  Modified on 8/30/2004 ####################
#		X:=Ord(denom(normal(x)),prime);
#		if X > 0 then 
#			RETURN(X/2);  # Point = Origin (mod prime)
#		fi;
		A:=Ord(numer(normal(3*x^2+2*a[2]*x+a[4]-a[1]*y)),prime)
			- Ord(denom(normal(3*x^2+2*a[2]*x+a[4]-a[1]*y)),prime);
		B:=Ord(numer(normal(2*y+a[1]*x+a[3])), prime)
			- Ord(denom(normal(2*y+a[1]*x+a[3])), prime);
		if A <=  0 or B <= 0 then 
			RETURN(max(0,Ord(denom(normal(x)),prime)/2));  # Point = Origin (mod prime)
		fi;
################################################################
		if member(prime, map2(op, 1, Bad_fibers), 'i') then
			if Bad_fibers[i][3] = 1 and  Bad_fibers[i][4] = 0 then  # Type I_b
				if Bad_fibers[i][5] = 1 then  # I_1
					RETURN(0);
				else
					B:=Ord(2*y+a[1]*x+a[3], prime);
					D:=Ord(Disc,prime);
					RETURN( min(B,D/2)*( min(B,D/2) - D ) / (2*D) );
				fi;
			elif Bad_fibers[i][3] = 2 then  # Type II or II*
				RETURN(0);
			elif Bad_fibers[i][3] = 4 then  # Type IV or IV*
				B:=Ord(2*y+a[1]*x+a[3], prime);
				RETURN(-B/3);
			else  # Type III, III* ou I_b*
				C:=Ord(3*x^4+b[2]*x^3+3*b[4]*x^2+3*b[6]*x+b[8], prime);
				RETURN(-C/8);
			fi;
		else
			RETURN(0);
		fi;

	elif prime = infinity then
		for i from 1 to 4 do a_inf[i]:=subs(t=1/t,a[i])*t^(i*(pg + 1)) od;
		for i from 1 to 4 do b_inf[2*i]:=subs(t=1/t,b[2*i])*t^(2*i*(pg + 1)) od;
		x_inf:=subs(t=1/t,point[1])*t^(2*(pg + 1));
		y_inf:=subs(t=1/t,point[2])*t^(3*(pg + 1));		
####################  Modified on 8/30/2004 ####################
#		X:=Ord(denom(normal(x_inf)),t);  
#		if X > 0 then 
#			RETURN(X/2);  # Point = Origin (mod prime)
#		fi;
		A:=Ord(numer(normal(3*x_inf^2+2*a_inf[2]*x_inf+a_inf[4]-a_inf[1]*y_inf)),t)
			- Ord(denom(normal(3*x_inf^2+2*a_inf[2]*x_inf+a_inf[4]-a_inf[1]*y_inf)),t);
		B:=Ord(numer(normal(2*y_inf+a_inf[1]*x_inf+a_inf[3])), t)
			- Ord(denom(normal(2*y_inf+a_inf[1]*x_inf+a_inf[3])), t);
		if A <= 0 or B <= 0 then 
			RETURN(max(0,Ord(denom(normal(x_inf)),t)/2));  # Point = Origin (mod prime)
		fi;
################################################################
		if member(prime, map2(op, 1, Bad_fibers), 'i') then
			if Bad_fibers[i][3] = 1 and  Bad_fibers[i][4] = 0 then  # Type I_b
				if Bad_fibers[i][5] = 1 then  # I_1
					RETURN(0);
				else
					B:=Ord(2*y_inf + a_inf[1]*x_inf + a_inf[3],t);
					D:=12*(pg + 1) - degree(Disc,t);
					RETURN( min(B,D/2)*( min(B,D/2) - D ) / (2*D) );
				fi;
			elif Bad_fibers[i][3] = 2 then  # Type II or II*
				RETURN(0);
			elif Bad_fibers[i][3] = 4 then  # Type IV or IV*
				B:=Ord(2*y_inf + a_inf[1]*x_inf + a_inf[3],t);
				RETURN(-B/3);
			else  # Type III, III* ou I_b*
				C:=Ord(3*(x_inf)^4+b_inf[2]*(x_inf)^3 + 3*b_inf[4]*(x_inf)^2
					 + 3*b_inf[6]*(x_inf) + b_inf[8],t);
				RETURN(-C/8);
			fi;
		else
			RETURN(0);
		fi;
	fi;
		
end:


Height:= proc(point)
	global a, b, c, Disc, j, pg, Bad_fibers;
	local places, denoms, i;
	description 
		`The canonical height of [x,y] with respect to the divisor (O).`;
	option remember;
	
	if point = `Origin` then 
		RETURN(0);
	fi;
	if Bad_fibers='not_yet' then
		Calculate_bad_fibers();
	fi;

	denoms:=map(L->[L[1],Repetition(L[1])],
				op(2,factors(denom(factor(point[1])))));
	places:= map(L->[L[1],L[2]],Bad_fibers);
	for i from 1 to nops(denoms) do
		if not member(denoms[i][1], map2(op,1,Bad_fibers)) then
			places:=Append(places,[denoms[i]])
		fi;
	od;
	if not member(infinity, map2(op,1,Bad_fibers)) then
		places:= Append(places,[[infinity,1]]);
	fi;
	( Add_all @ map2 ) 
		((pt,pl) -> Local_height(pt,pl[1])*pl[2], point , places)
	+ (pg + 1) ;
end:


Height_pairing:= proc(p1, p2)
	global a, b, c, Disc, j, pg, Bad_fibers;
	description 
		`The canonical height pairing of the two points, \
		normalized so that Height(P) equals Height_pairing(P,P).`;
	option remember;
		
	if p1 = p2 then Height(p1) 
	else
		1/2 * ( Height(p1 &++ p2) - Height(p1) - Height(p2) );
	fi;
end:

Height_matrix:= proc()
	global a, b, c, Disc, j, pg, Bad_fibers;
	local M, m, n;
	description `The matrix of height pairings for the given set of points.`;
	
	M:= array( symmetric, 1..nargs,1..nargs ):
	for m from 1 to nargs do
		for n from 1 to m do
			M[m,n]:=Height_pairing(args[m],args[n]);
		od;
	od;
	op(M): # I haven't quite understood this, but it works.
end:

Height_regulator:= proc()
	global a, b, c, Disc, j, pg, Bad_fibers;
	description `The height regulator for the given set of points. \
		(N.B. This differs from the regulator in the Birch-Swinnerton-Dyer \
		conjecture by a power of 2.)`;
	
	linalg[det](Height_matrix(args));
end:

################################################################
#                    Division Polynomials                      #
################################################################

Division[3]:=3*X^4+(a[1]^2+4*a[2])*X^3+3*(a[1]*a[3]+2*a[4])*X^2 
+3*(a[3]^2+4*a[6])*X+1/4*(a[1]^2+4*a[2])*(a[3]^2+4*a[6]) 
-1/4*(a[1]*a[3]+2*a[4])^2:

# Division[2n] differs usual psi_2n by a factor of 2*y+a[1]*x+a[3].
Division[4]:=2*X^6+(a[1]^2+4*a[2])*X^5+5*(a[1]*a[3]+2*a[4])*X^4
+10*(a[3]^2+4*a[6])*X^3 +10*(1/4*(a[1]^2+4*a[2])*(a[3]^2+4*a[6])
-1/4*(a[1]*a[3]+2*a[4])^2)*X^2
+((a[1]^2+4*a[2])*(1/4*(a[1]^2+4*a[2])*(a[3]^2+4*a[6])
-1/4*(a[1]*a[3]+2*a[4])^2)-(a[1]*a[3]+2*a[4])*(a[3]^2+4*a[6]))*X
+(a[1]*a[3]+2*a[4])*(1/4*(a[1]^2+4*a[2])*(a[3]^2+4*a[6])
-1/4*(a[1]*a[3]+2*a[4])^2)-(a[3]^2+4*a[6])^2:



################################################################
#                    Utility functions                         #
################################################################

Ord:= proc ()  
	local exp, prime, point, n;
	
	if nargs = 2 then 
		# ord(poly, prime) = how many times `prime` divides `exp`.
		# !!! does not work well over alg num field !!!

		exp:=normal(args[1]);
		prime:=args[2];
		
		if  type(exp, polynom) then
			if exp = 0 then 
				RETURN(infinity)
			else 
				(Add_all @ map) (
					(x,y) -> if y = x[1] then x[2] else 0 fi,
					op(2,factors(exp)),
					args[2]
				);
			fi;
	 	else 
			(Add_all @ map) (
				(x,y) -> if y = x[1] then x[2] else 0 fi,
				op(2,factors(numer(normal(exp)))),
				prime
			)
			- (Add_all @ map) (
				(x,y) -> if y = x[1] then x[2] else 0 fi,
				op(2,factors(denom(normal(exp)))),
				prime
			);		
		fi;
		
	elif nargs = 1 and type(args[1],list) and nops(args[1]) = 2 then
		# calculate the order of a point.
		
		point:=args[1];
		if Height(point) <> 0 then 
			RETURN(infinity);
		else
			n:=1;
			while point <> `Origin` do
				point:= point &++ args[1];
				n:= n + 1;
			od;
			RETURN(n);
		fi;
		
	elif nargs = 1 and args[1] = `Origin` then
		RETURN(1);
	
	else
		ERROR(`wrong type of arguments in function Ord.`);
	
	fi;
end:

Append:= proc (list1, list2)
	description `Append list2 at the end of list1.`;
	[op(list1),op(list2)]
end:

Add_all:= proc (list)
	description `Take the sum of all elements in the list.`;
	( `+` @ op ) (list);
end:

Multiply_all:= proc (list)
	description `Take the product of all elements in the list.`;
	( `*` @ op ) (list);
end:
	
		
macro(
a=a,
b=b,
c=c,
j=j,
Disc=Disc,
pg=pg,
t=t,
curve=curve):

#
# This definition must come after closing the macro.
#
Show:=proc(key,opt)
	global a, b, c, Disc, j, pg, Bad_fibers;

	if  key='curve' then 
	    if nargs<2 or opt='original' or opt='long' then
		     RETURN(y^2 + `«a»`[1] * x * y + `«a»`[3] * y 
			     = sort(x^3 + `«a»`[2] * x^2 + `«a»`[4] * x + `«a»`[6],x));
		elif opt='homogeneous' or opt='homg' or opt='homo' or opt='homog' then 
		     RETURN(y^2*z + `«a»`[1] * x * y * z + `«a»`[3] * y *z^2
			     = sort(x^3 + `«a»`[2] * x^2 * z + `«a»`[4] * x *z^2 + `«a»`[6] *z^3, x));
		elif nargs=2 and ( opt='homogeneous' or opt='homg' or opt='homo' ) then 
		     RETURN(y^2*z + `«a»`[1] * x * y * z + `«a»`[3] * y *z^2
			     =  x^3 + `«a»`[2] * x^2 * z + `«a»`[4] * x *z^2 + `«a»`[6] *z^3, x);
		elif nargs=2 and ( opt='short' or opt='simplified' ) then
		     RETURN(y^2 
			     = sort(x^3 - 27* `«c»`[4] * x - 54* `«c»`[6],x));
		elif nargs>2 and ( opt='short' or opt='simplified' or opt='homogeneous') then
		     # à amérioler
		     RETURN(y^2 
			     = x^3 - 27* `«c»`[4] * x * z^2 - 54* `«c»`[6] * z^3);
		fi:
	elif key='Disc' then  RETURN(`«Disc»`): 
	elif key='j' or  key='jay' then  RETURN(`«j»`): 
	elif key='pg' then RETURN(`«pg»`):
	elif key='`bad fibers`' or key='bad_fibers' then Show_bad_fibers():
	elif key='all' or key='All' then Show_data():
	elif key='a' then RETURN([`«a»`[1],`«a»`[2],`«a»`[3],`«a»`[4],`«a»`[6]]):
	elif key='b' then RETURN([`«b»`[2],`«b»`[4],`«b»`[6],`«b»`[8]]):
	elif key='c' then RETURN([`«c»`[4],`«c»`[6]]):	
	fi; 
end:

################################################################
#                    Added in April 2004                       #
################################################################
	
Homogenize:=proc(exp)
	description "Convert a polynomial in x, y or in X, Y to a homogeneous polynomial";
	if degree(exp, {x,y}) > 0 then
		map(factor,collect(subs({x=x/z,y=y/z},exp)*z^degree(exp,{x,y}),{x,y,z}));
	elif degree(exp, {X,Y}) > 0 then
		map(factor,collect(subs({X=X/Z,Y=Y/Z},exp)*Z^degree(exp,{x,y}),{X,Y,Z}));
	else exp;
	fi:
end:

## `conversion` is the result of CtoW
Trcw:=proc(curve,trwc,trcw,point)
	factor(subs({x=point[1],y=point[2],z=point[3]},
	[trcw[1]/trcw[3],trcw[2]/trcw[3]]));
end:
	
################################################################
#                    Added in June 2008                        #
################################################################
mapfactor:=proc(expression,variables) 
          local nonconst_part, const_part;
		nonconst_part:=simplify(expression
					- subs(map(x->x=0,variables),simplify(expression)));
		nonconst_part:=map(factor, collect(nonconst_part,
 				variables, 'distributed'));
 		const_part:= factor(subs(map(x->x=0,variables),simplify(expression)));
		RETURN(nonconst_part+const_part);
end: