# $Source: /u/maple/research/lib/algcurves/src/RCS/Weierstrassform,v $
# $Notify: mvanhoei@daisy.uwaterloo.ca $

macro(
	RAD={radical,nonreal},
	Weierstrassform=`algcurves/Weierstrassform`,
	find_a_point=`algcurves/find_a_point`,
	find_points=`algcurves/find_points`,
	degree_ext=`algcurves/degree_ext`,
	g_conversion1=`algcurves/g_conversion1`,
	g_conversion2=`algcurves/g_conversion2`,
	g_factors=`algcurves/g_factors`,
	g_ext=`algcurves/g_ext`,
	g_ext_r=`algcurves/g_ext_r`
):

Weierstrassform:=proc(F,x::name,y::name,x0::name,y0::name)
	global `algcurves/x0`,`algcurves/y0`;
	local f,p,g,d,P,D;
	options remember, 
`Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
	# Make sure that the coefficients are evala normalized:
	Normalizer:=`if`(has([args],RootOf),`evala/Normal`,normal);
	f:=collect(primpart(F,{x,y}),{x,y},'distributed',Normalizer);
if indets([args],RAD)<>{} then
	f:=[args];
	p := radfield(indets(f,RAD));
	return eval(subs(p[2],procname(op(eval(subs(p[1],f)))))) 
elif not has(f,x) or not has (f,y) then
	error "%1 must have both variables %2,%3",F,x,y
elif has(evala(Content(f,x)),y) or has(evala(Content(f,y)),x) then
	error "%1 is reducible",F
elif nargs<6 or not type(args[6],list) then
	p,g:=find_a_point(f,x,y);
    if g<>2 then
	d:=degree_ext(p,[args]);
	if d>1 then
		# Try to find a better point.
		P:=find_a_point(f,y,x)[1];
		D:=degree_ext(P,[args]);
		if D<d or (D=d and length(P)<length(p)) then
			p:=[P[2],P[1],P[3]]
		fi;
		if has(p,RootOf) then
			Normalizer:=`evala/Normal`
		fi
	fi;
	userinfo(2,'algcurves',`Using point`,p)
    fi
else
	g:=`algcurves/genus`(F,x,y);
	p:=NULL
fi;
if g<1 then
	error "Not (hyper)-elliptic, genus is",g
elif g=1 then
	# The reason for using this global `algcurves/x0` is that the
	# options remember in `algcurves/genus1` will be more effective
	# this way.
	subs(`algcurves/x0`=x0,`algcurves/y0`=y0,`algcurves/genus1`(
	 f,x,y,`algcurves/x0`,`algcurves/y0`,p,args[6..nargs]))
else
	`algcurves/genus2`(f,x,y,x0,y0,g,p,args[6..nargs])
fi
end:


# Find a regular point on the curve
# For use only in the procedure Weierstrassform, for taking
# a regular point and checking g=1.
find_a_point:=proc(F,x,y)
	local vp,i,ext,f,d,m,g;
	global g_conversion1,g_conversion2;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
	ext:=g_ext(F);
	f:=subs(g_conversion1,F);
	vp:=find_points(f,x,y,ext,[`search points on`,{1,2,3}]);
	d:=degree(f,{x,y});
	g:=(d-1)*(d-2)/2-add(i[3]*degree_ext(i,f),i=vp);
	# Add a regular point, just in case we found no regular points yet.
	i:=0;
	while testeq(discrim(subs(x=i,F),y)) do i:=i+1 od;
	vp:=vp union {seq([[i,RootOf(d[1],y),1],1,0],
	 d=g_factors(subs(x=i,f),y,g_ext_r(F)))};
	d:={};
	m:=min(seq(`if`(i[2]=1,degree_ext(i,f),infinity),i=vp));
	for i in vp do if degree_ext(i,f)=m and i[2]=1 then
		d:=d union {i[1]}
	fi od;
	# Take the smallest one.
	m:=min(seq(length(i),i=d));
	for i in d do if length(i)=m then
		return subs(g_conversion2,i),g
	fi od;
	# The following can't happen, d can not be empty because
	# there should at least be one regular point in there (we
	# placed it there in the beginning of this procedure.
	error "This is not supposed to happen"
end:

#savelib ('Weierstrassform','find_a_point'):