# $Source: /u/maple/research/lib/algcurves/src/RCS/singularities,v $
# $Notify: hoeij@sci.kun.nl $

macro(
	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`. These files will
# be read in the singularities procedure

macro(
	homogeneous=`algcurves/homogeneous`
):
macro(
	local_intbasis23=`algcurves/ib23`,
	double_factors=`algcurves/db_factors`,
	local_intbasis=`algcurves/local_ib`,
	ext_to_coeffs=`algcurves/e_to_coeff`,
	g_gcdex=`algcurves/gcdex`
):

macro(
	v_ext_m=`algcurves/v_ext_m`,
	lift_exp=`algcurves/lift_exp`,
	lift_exp_m1=`algcurves/lift_exp_m1`,
	truncate_subs=`algcurves/truncate_subs`,
	monic=`algcurves/monic`,
	`puiseux/technical_answer`=`algcurves/puiseux_te`,
	`integral_basis/bound`=`algcurves/ib_bound`,
	Newtonpolygon=`algcurves/Newtonpolygon`
):

macro(
	g_conversion1=`algcurves/g_conversion1`,
	g_conversion2=`algcurves/g_conversion2`,
	g_normal=`algcurves/g_normal`,
	g_expand=`algcurves/g_expand`,
	normal_tcoeff=`algcurves/normal_tcoeff`,
	g_evala=`algcurves/g_evala`,
	g_evala_rem=`algcurves/g_evala_rem`,
	g_zero_of=`algcurves/g_zero_of`,
	g_factors=`algcurves/g_factors`,
	rootof=`algcurves/rootof`,
	g_ext=`algcurves/g_ext`,
	g_ext_r=`algcurves/g_ext_r`,
	truncate=`algcurves/truncate`
):

# 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,k,mult,pl_transl,res;
	options remember, 
	    `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	n:=degree(f,{x,y});
if nargs=4 or (nargs=5 and has(args[5],`search points on`)) then
	j:=discrim(f,y);
	# First a quick testeq, then if necessary a thorough test g_normal=0
	if testeq(subs(g_conversion2,j)) and g_normal(j)=0 then
		ERROR(`Input should be squarefree`)
	fi;
	{seq(op(procname(args[1..4],i,args[5..nargs])),i=
	 [op(double_factors(j,x,ext)),infinity])}
elif where=infinity then
	k:=-1;
	while k=-1 or degree(subs(z=0,ff),x)<n do
		# Making sure that we may assume y=1
		k:=k+1;
		ff:=homogeneous(subs(y=y+k*x,f),x,y,z,polynom)
	od;
	res:={seq([[i[1][2],1+k*i[1][2],0],op(2..nops(i),i)],
	 i=procname(subs({z=x,x=y},subs(y=1,ff)),x,y,ext,[x,4]))};
	if k<>0 and not member([0,0],{seq([i[1][2],i[1][3]],i=res)}) then
		res union {[[1,0,0],1,0]}
	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])) then
		i:=local_intbasis23(f,x,y,ext,p,op(where),z);
		if not type(i,list) then
			# i should be of degree 1 in y
			subs(g_conversion1,
			 {[[p,g_normal(-coeff(i,y,0)/coeff(i,y,1)),1],2,1]})
		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);
		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;
			for j in pl_transl do if subs(x=0,j[1])=i then
			 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]
		od;
		{res}
	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.`;
if nargs<3 then
	ERROR(`wrong number of arguments`)
elif indets([args],radical)<>{} then
	F:=[args];
	F:=procname(op(subs(readlib(`radnormal/radfield`)(F,'v'),F)));
	RETURN(subs(eval(v),F))
elif not has(f,{x,y}) then
	ERROR(`Input should be a polynomial in`,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;

readlib(`algcurves/g_expand`):
readlib(`algcurves/puiseux`):
readlib(`algcurves/integral_basis`):
readlib(`algcurves/homogeneous`):
	ext:=g_ext([args]);
	F:=g_expand(subs(g_conversion1,f),ext);
	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','`algcurves/singularities.m`'):