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

macro(
	genus2=`algcurves/genus2`,
	index2subfield=`algcurves/index2subfield`
):

# Examples of curves with genus 2:
# f:=x^5+y^5+x^4*y^3;
# f:=y^7+(x-1)^3*x^4+(x-1)^3*x^2*y^2;
# Weierstrassform(f,x,y,x0,y0);

macro(
	EXT = {nonreal, radical, algext},
	E1 = `curve is a square`
):
# quotes " instead of ` don't work here with E1 in Release 5,
# cf. test file f62
# Cause:
# h:=traperror( ERROR( "hallo daar" ));
# evalb(h="hallo daar");  -> false in R5, true in next version.

genus2:=proc(f,x::name,y::name,x0::name,y0::name,g::posint,pt::list)
	local v,i,a,b,d,nd,dd,ni,f0,x0v,y0v,yv,id,z;
	global skip_dx;
	option
`Copyright (c) 1999 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
	ni:=member(`no inverse`,{args});
	if degree(f,y)=2 then
		# x = x0
		d:=Normalizer(discrim(f,y)/4);
		if d=0 then
			ERROR( E1 )
		fi;
		nd:=denom(d);
		d:=sqrfree(Normalizer(d*nd^2),x);
		a:=[sqrt(d[1],symbolic),sqrt(-d[1],symbolic)];
		b:=indets(d[1],EXT);
		if indets(a[1],EXT)=b then
			a,b:=a[1],1
		elif indets(a[2],EXT)=b then
			a,b:=a[2],-1
		else
			a,b:=1,d[1]
		fi;
		dd:=a*mul(i[1]^iquo(i[2],2),i=d[2])/nd;
		d:=b*mul(`if`(irem(i[2],2)=1,i[1],1),i=d[2]);
		# Now discrim = d * dd^2,  d=squarefree.
		b,a:=seq(coeff(f,y,i),i=1..2);
		v:=[y0^2-subs(x=x0,d),x,(a*y+b/2)/dd,op(`if`(ni,
		 [],[x0,subs(x=x0,(-b/2+dd*y0)/a)]))];
		if degree(d,x)=2*g+1 then
			RETURN( v )
		elif degree(d,x)<>2*g+2 then
			ERROR( "should not happen" )
		fi;
		a:=roots(d,x,indets(f,RootOf));
		if a=[] then
			RETURN( v )
		fi;
		# Pick smallest root:
		a:=sort(a,(S,T) -> evalb(length(S)<length(T)))[1][1];
		b:=y0^2-Normalizer(x0^(2*g+2)*subs(x=1/x0+a,d));
		if g>1 and length(b)>3/2*length(v[1])+5 then
			userinfo(3,'algcurves',`Skipping reducing degree by 1`);
			RETURN( v )
		fi;
		userinfo(3,'algcurves',`Reducing degree by 1`);
		RETURN( map(Normalizer,[b,1/(x-a),v[3]/(x-a)^(g+1),
		 op(subs(x0=1/x0+a,y0=y0/x0^(g+1),v[4..nops(v)]))]) )
	elif degree(f,x)=2 then
		v:=procname(f,y,x,x0,y0,g,[],args[8..nargs]);
		RETURN( [v[1],v[2],v[3],op(`if`(ni,
		 [],[v[5],v[4]]))] )
	elif degree(f,{x,y})-ldegree(f,{x,y})=2 then
		v:=procname(subs(x=1,`algcurves/homogeneous`(f,x,y,z))
		 ,y,z,x0,y0,g,[],args[8..nargs]);
		RETURN( map(Normalizer,[v[1],op(subs(y=y/x,z=1/x,[v[2],v[3]]))
		 ,op(`if`(ni,[],[1/v[5],v[4]/v[5]]))]) )
	fi;
	if member(E1,{args}) then
		# x,y have interchanged.
		b:=args[nargs]
	elif degree(f,x)>2 and degree(f,y)>2 then
		b:=`algcurves/differentials`(f,x,y,skip_dx);
		if nops(b)<>g then
			ERROR( "curve must be reducible" )
		elif g>2 then
			b:=index2subfield(f,x,y,g,b,map(Normalizer,pt),args[8..nargs])
		fi
	else
		ERROR( "should not happen" )
	fi;
	v:=Normalizer(b[1]/b[2]);
	v:=map(primpart,[numer(v),denom(v)]);
	if (length(v[1])<length(v[2]) and has(v[2],y)) or not has(v[1],y) then
		v:=[v[2],v[1]]
	fi;
	x0v:=v[1]/v[2];
	# x0 * v[1] - v[2] = 0  -> leaves 2 values for y.
	# v[1]=Equation for y, v[2]=eqn for x:
# v:=[seq(numer(evala(Normal(
# evala(Resultant(x0*v[2]-v[1],f,i))/evala(Resultant(v[2],f,i)),expanded
# ))),i=[x,`if`(ni,NULL,y)])];
# Doesn't work, see:
# f:=
# 20-28*x*y^3+16*x+7*y^6+28*x^2-2*y^2-2*x^4+2*x^3-y^5+y^3*x^2-8*y*x^3+4*y^4*x^2+
# 2*y^2*x-y*x^5+x^3*y^3-20*y*x^2+5*y^2*x^4-8*y^3;
# Weierstrassform(f,x,y,x0,y0); -> error, "degree is not 2"
#
# The second one Resultant(v[2],f,i) should have been Resultant(v[1],f,i), but I
# think the following is faster:
	v:=[seq(evala(Primpart(
	 evala(Primpart(evala(Resultant(x0*v[2]-v[1],f,i)),x0)),{x,y}))
	,i=[x,`if`(ni,NULL,y)])];
	if degree(v[1],y)<>2 or ((not ni) and degree(v[2],x)<>2) then
		ERROR( "degree is not 2" )
	fi;
	a:=traperror(procname(v[1],x0,y,x0,y0,g,pt));
	if a=lasterror then
		if a=E1 and not member(a,{args}) then
			v:=procname(f,y,x,x0,y0,
			 g,[pt[2],pt[1],pt[3]],args[8..nargs],a,b);
			RETURN( [v[1],v[2],v[3],op(`if`(ni,
			 [],[v[5],v[4]]))] )
		else
			ERROR( a )
		fi
	fi;
	f0,x0v,y0v,yv:=a[1],Normalizer(subs(x0=x0v,a[2])),
	 `algcurves/FFDIV`(f,x,y,Normalizer(subs(x0=x0v,a[3]))),a[5];
	if ni then
		RETURN( [f0,x0v,y0v] )
	fi;
	v:=numer(Normalizer(subs(x0=a[4],v[2])));
	b,a:=seq(coeff(v,x,i),i=1..2);
	d:=Normalizer( discrim(v,x)/(y0^2-f0) ); # must be a square
	if d=0 then
		RETURN( [f0,x0v,y0v,-b/2/a,yv] )
	fi;
	d:=sqrfree(d,x);
	for i in d[2] do if type(i[2],odd) then
		ERROR( "not a square" )
	fi od;
	id:=indets({args},EXT);
	d:=mul(i[1]^iquo(i[2],2),i=d[2])*RootOf(evala(Factors(
	 x^2-d[1],id))[2][1][1],x);
	# xv = (-b +/- d*y0)/2a
	# y0 = (2a*x+b)/(+/- d), so one of these must be zero:
	v:=[y0 - (2*a*x+b)/d, y0 + (2*a*x+b)/d];
	v:=map(`algcurves/homogeneous`,map(numer@Normalizer,v),x0,y0,z);
	i:=lcm(denom(x0v),denom(y0v));
	nd:=x0=Normalizer(x0v*i), y0=Normalizer(y0v*i), z=i;
	for i from 0 to 100 do
		userinfo(5,'algcurves',`zero testing...`);
		dd:=evala(Factors(subs(x=i,f),id))[2];
		if not has(dd,y) then next fi;
		dd:=map(evala@Normal,subs(nd,x=i,y=
		 RootOf(dd[1][1],y),v));
		if dd=[0,0] then next fi;
		if dd[1]=0 then
			RETURN( [f0,x0v,y0v,(-b + d*y0)/2/a,yv] )
		elif dd[2]=0 then
			RETURN( [f0,x0v,y0v,(-b - d*y0)/2/a,yv] )
		else
			break
		fi
	od;
	ERROR( "should not reach this point" )
end:


# f = an irreducible polynomial in x,y
# g = the genus. g>2.
# B = algcurves[differentials](f,x,y,skip_dx)
# pt = a list [x,y,z] giving a regular point on the curve
#      in homogeneous coordinates.
index2subfield:=proc(f,x,y,g,B,pt,var)
	local z,su,a,i,j,k;
	option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
if type(B,list) then
	if not `algcurves/is_hyperelliptic`(f,x,y) then
		error "Not (hyper)-elliptic, genus is %1", g 
	fi;
	z:={seq(a[i],i=1..g)};
	k:=add(a[i]*B[i],i=1..g);
	j:=g;
	while nops(z)>2 do
		k:=procname(f,x,y,j,k,pt,z);
		z:=indets(k) intersect z;
		j:=j+nops(z)-2
	od;
	[subs(z[1]=0,z[2]=1,k),subs(z[1]=1,z[2]=0,k)]
elif pt[3]=0 then
	if not type(pt,list) or nops(pt)<>3 or pt=[0,0,0] then
		ERROR( "wrong arguments" )
	fi;
	k:=subs(x=1,map(`algcurves/homogeneous`,[f,B],x,y,z));
	k:=procname(k[1],y,z,g,k[2],[pt[2],pt[3],pt[1]],var);
	subs(z=1,`algcurves/homogeneous`(k,y,z,x))
else
	su:=x=pt[1]/pt[3],y=pt[2]/pt[3];
	if Normalizer(subs(su,f))<>0 then
		ERROR( "not a point on the curve" )
	elif Normalizer(subs(su,diff(f,y)))=0 then
		if Normalizer(subs(su,diff(f,x)))=0 then
			ERROR( "given point is a singularity" )
		else
			RETURN( procname(f,y,x,g,B,[pt[2],pt[1],pt[3]],var) )
		fi
	fi;
	# Can be done more efficiently by Puiseux expansions and by using
	# linear equations from other points as well.
	j:=subs(su[1],f);
	j:=evala(Gcd(j,diff(j,y)));
	subs(`solve/linear`({seq(subs(x=pt[1]/pt[3],Normalizer(
	 subs(y=pt[2]/pt[3],`if`(i=0,B,`algcurves/FFDIV`(f,x,y,
	  normal(subs(RootOf(f,y)=y,j^(2*i-1)*diff(subs(y=RootOf(f,y),B),x$i))))
	)))),i=0..g-3)}, var),B)
fi
end:

#savelib('genus2','index2subfield','`algcurves/genus2.m`'):