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

macro(
	RAD={radical,nonreal},
	singularities=`algcurves/singularities`,
	find_points=`algcurves/find_points`,
	degree_ext=`algcurves/degree_ext`
):

# This code uses the procedures in `algcurves/g_expand`,
# `algcurves/puiseux` and `algcurves/integral_basis`.

macro(
	homogeneous=`algcurves/homogeneous`,
	local_intbasis23=`algcurves/ib23`,
	double_factors=`algcurves/db_factors`,
	lift_exp=`algcurves/lift_exp`,
	`puiseux/technical_answer`=`algcurves/puiseux_te`,
	g_conversion1=`algcurves/g_conversion1`,
	g_conversion2=`algcurves/g_conversion2`,
	g_normal=`algcurves/g_normal`,
	g_expand=`algcurves/g_expand`,
	g_zero_of=`algcurves/g_zero_of`,
	g_ext=`algcurves/g_ext`
):

# This is used as a global variable for making the curve homogeneous. The
# reason for using a global variable here is that the options remembers in
# various places in this package will be more useful this way.
macro(
	z=`algcurves/z`
):

# Input: f a polynomial (in expanded form) in x and y.
# rootof syntax
# Output: the singularities, plus some other info that is used
# in odd_singularity_on_C
# Furthermore: the multiplicity and the delta invariant of these points
find_points:=proc(f,x,y,ext,where)
	global g_conversion1,z;
	local ff,n,p,delta,extp,f_translated,i,j,mult,pl_transl,res,r;
	options remember, 
`Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
	n:=degree(f,{x,y});
if nargs=4 or (nargs=5 and has(args[5],`search points on`)) then
	{seq(op(procname(args[1..4],i,args[5..nargs])),i=
	 [op(double_factors(f,x,y,ext)),infinity])}
elif where=infinity then
	ff:=homogeneous(f,x,y,z,polynom);
	res:={seq(subsop(1=[i[1][2],1,0],i),
	 i=procname(subs({z=x,x=y},subs(y=1,ff)),x,y,ext,[x,4]))};
	p:=n-degree(subs(z=0,ff),x);
	if p>0 then
		mult:=0;
		delta:=0;
		r:=0;
	    if p>1 then
		pl_transl:=[op(lift_exp([0,1,z,ext,p,1,0,[]],subs(x=1,ff),z,y))];
		for j in pl_transl do
		 r:=r+degree_ext(j,ext);
		 mult:=mult+degree_ext(j,ext)*min(degree(j[3],z),
		  max(ldegree(j[1],z),1));
		 delta:=delta+1/2*degree_ext(j,ext)*
		  (degree(j[3],z)*j[7]-degree(j[3],z)+1)
		od;
		if not (type(mult,integer) and type(delta,integer))
			then error "found wrong values"
		fi;
	    fi;
		res union {[[1,0,0],max(1,mult),delta,r]}
	else
		res
	fi
elif where[2]<2 then
	{}
else
	p:=g_zero_of(where[1],x,'extp');
	extp:=eval(extp);
	if member(where[2],{2,3}) and not
	# If the following is the case, then use the other method because
	# the other method (although slower) will be efficient enough
	# (because degree_ext is small) and with some luck it may give
	# an odd degree place.
	(has([args],`search points on`) and member(degree_ext(p,ext),
	 args[6][2])) and g_normal(subs(x=p,lcoeff(f,y)))<>0 then
		i:=local_intbasis23(f,x,y,ext,p,op(where),z);
		if nops([i])=2 then
			# i should be of degree 1 in y
			subs(g_conversion1,
			 {[[p,g_normal(-coeff(i[1],y,0)/coeff(i[1],y,1)),1]
			,2,1, `if`(i[2],1,2)]})
		else
			# perhaps this line is useful in odd_singularity_on_C, if
			# where[2]=3 we have a place of algebraic odd-degree 
			{[[p,0,1],0,0,where[2]]}
		fi
	else
		res:=NULL;
		extp:=[extp,op(ext)];
		f_translated:=g_expand(subs(x=x+p,f),extp);
		i:=subs(x=0,f_translated);
		j:=degree(i,y);
	    if j=0 then
		{}
	    elif j=1 then
		{[[p,Normalizer(-coeff(i,y,0)/coeff(i,y,1)),1],1,0]}
		# The case j=1 had to be handled seperately; it can't be
		# handled by the method below because the puiseux code won't
		# do any lifting steps when j=1.
	    else
		pl_transl:=[op(`puiseux/technical_answer`(f_translated,x,y
		 ,0,extp))];
		if irem(degree_ext(p,ext),2)=1 then
			# add info that could be useful in odd_singularity_on_C
			for i in pl_transl do
			if irem(i[6],2)=1 then
				res:=[[p,0,1],0,0,3];
				break
			fi od
		fi;
		ff:={seq(subs(x=0,i[1]),i=pl_transl)}; # the y values
		for i in ff do # i is y value
			mult:=0;
			delta:=0;
			r:=0;
			for j in pl_transl do if subs(x=0,j[1])=i then
			 r:=r+degree_ext(j,[i,op(extp)]);
			 mult:=mult+degree_ext(j,[i,op(extp)])*min(
			 # Minimum of the valuation of x and y in this place
			  degree(j[3],x),
			  max(ldegree(j[1]-subs(x=0,j[1]),x),1)
			 );
			 delta:=delta+1/2*degree_ext(j,[i,op(extp)])*
			  (degree(j[3],x)*j[7]-degree(j[3],x)+1)
			fi od;
			if not (type(mult,integer) and type(delta,integer))
				then error "found wrong values"
			fi;
			res:=res,[[p,i,1],mult,delta,r]
		od;
		{res}
	    fi
	fi
fi
end:

# Input: f is a polynomial in x and y
# Output is a list of the following lists:
# [ [x,y,z], multiplicity, delta invariant ]
singularities:=proc(f,x,y)
	global g_conversion1,g_conversion2;
	local F,v,i,res,ext;
	options remember, 
`Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
if nargs<3 then
	error "wrong number of arguments"
elif indets([args],RAD)<>{} then
	F:=[args];
	v := radfield(indets(F,RAD));
	return eval(subs(v[2],procname(op(eval(subs(v[1],F)))))) 
elif not has(f,{x,y}) then
	error "Input should be a polynomial in %1 %2",x,y
elif has(evala(Content(f,y)),x) then
	# If f has a component which is independent of y, this component
	# will be ignored by `algcurves/find_points`, but I want to
	# include it in the computation. This is a fix for this case
	v:=procname(subs(x=x+y,f),x,y,args[4..nargs]);
	return {seq([[i[1][1]+i[1][2],i[1][2],i[1][3]],i[2],i[3]],i=v)}
fi;

	ext:=g_ext([args]);
	F:=subs(g_conversion1,
	 collect(primpart(f,{x,y}),{x,y},'distributed',Normalizer));
	v:=find_points(F,x,y,ext);
	res:=NULL;
	for i in v do if i[2]>1 then
		res:=res,i
	fi od;
	subs(g_conversion2,{res})
end:

#savelib ('find_points','singularities'):