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

macro(	solve=readlib(`solve/linear`)	);

macro(	iss94=`algcurves/iss94`,
	function_with_one_pole=`algcurves/f_with_1_p`
):

macro(
	homogeneous=`algcurves/homogeneous`
):
macro(
	singularities=`algcurves/singularities`,
	find_points=`algcurves/find_points`,
	degree_ext=`algcurves/degree_ext`
):
macro(
	integral_basis=`algcurves/integral_basis`,
	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(
	puiseux=`algcurves/puiseux`,
	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`
):
macro(	ratpar=`algcurves/ratpar`,
	odd_point_on_a_old=`algcurves/oddp_a`,
	inverse_of_g=`algcurves/inv_g`,
	genus1_charpol=`algcurves/g1_charpol`,

	odd_root=`algcurves/odd_root`,
	regpoint_from_sing=`algcurves/rp_from_s`,
	odd_regpoint_C=`algcurves/oddrp_C`,
	odd_point_on_C2=`algcurves/oddp_C`,
	rat_point_on_C2=`algcurves/ratp_C`,
	search_rat_param=`algcurves/s_param`,
	odd_singularity_on_C=`algcurves/odds_C`,
	parametrize_cube=`algcurves/param_cube`,
	parametrize_conic=`algcurves/param_conic`,

	express_in_p=`algcurves/expr_in_p`,
	express_x_in_p=`algcurves/expr_x_in_p`,
	compute_x_old=`algcurves/comp_x`,
	frac_integral_a=`algcurves/frac_int_a`,
	FFDIV=`algcurves/FFDIV`,
	L_inf=`algcurves/L_inf`
):

macro(
	z=`algcurves/z`
):

# f is a polynomial in x and y
# f must be irreducible over Qbar
# If f is reducible this procedure does not work properly.
iss94:=proc(f,x,y,t,point)
	global z;
	local d,fu;
	options remember, 
	    `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	d:=degree(f,{x,y});
if coeff(f,y,d)=0 then
	if coeff(f,x,d)<>0 then
		# switch x and y
		d:=procname(f,y,x,t,[point[2],point[1],point[3]]);
		RETURN([d[2],d[1]])
	elif coeff(coeff(f,x,0),y,0)<>0 then
		# switch z and y
		d:=procname(subs(y=1,z=y,homogeneous(f,x,y,z,
		 polynom)),x,y,t,[point[1],point[3],point[2]]);
		RETURN(map(g_normal,[d[1]/d[2],1/d[2]]))
	else
		d:=procname(collect(subs(x=x+y,f),{x,y},'distributed',g_normal)
		 ,x,y,t,[point[1]-point[2],point[2],point[3]]);
		RETURN([g_normal(d[1]+d[2]),d[2]])
	fi
fi;
	# coeff(f,y,d)<>0, so [0,1,0] is not a point on the curve anymore.
	# Now we compute a function that has a pole in point. If the
	# point is finite, this function has no other poles in finite
	# otherwise it has no other poles in infinity.
	# It may still have other poles, they will be removed by
	# the procedure function_with_one_pole
	if point[3]=0 then
		fu:=[(homogeneous(evala(Quo(subs({t=0,x=1}
		 ,homogeneous(f,x,y,t,polynom)),y-point[2]/point[1],y))
		 ,y,t,x,polynom)),infinity];
	else
		fu:=[evala(Quo(evala(Expand(subs(x=point[1]/point[3],f)))
		 ,y-point[2]/point[3],y))/(x-point[1]/point[3]),'finite']
	fi;
	fu:=function_with_one_pole(f,x,y,fu);
	userinfo(2,'algcurves',`Computed a generator of the function field`);
	express_in_p(f,x,y,t,numer(fu),denom(fu))
end:


# The following procedure computes a parameter. A parameter is a function
# that has only 1 pole, with multiplicity 1. Then k(p)=k(x,y), p generates
# the function field. The pole will be in point. In fact, as the syntax
# is now, point is not a point, but a function with a pole in a certain point.
#
# point=[function,'finite'] or [function,infinity]. This function must
# have one pole in finite or infinity. This function is the start for our
# search for a parameter. We add an ansatz to this point. Since this function
# is already good in a part of the plane (either 'finite' or infinity)
# we will change no poles in that part of the plane, and try to remove the
# poles in the other part.
#
# The point [0,1,0] must not be a point on the curve. The reason is that
# we have divided the projective plane in 2 affine parts, a line and
# a plane. Doing this, there remains 1 point, and therefore that point must
# not be on the curve.
#
# Note that the way we divide the plane in parts is not really relevant
# for the method. The fact that one part is only a line does not give much
# asymmetry in the algorithm. The important thing about these parts is that
# they are affine, and not projective. We need this to apply integral basis.
#
# We will compute an integral basis for the finite part of the plane, and a
# local integral basis for the line infinity. Using these integral basis
# we can find linear equations stating that an ansatz, or an ansatz + a
# given function is integral (i.e. no poles) in a part of the plane (the
# part 'finite' and infinity.
function_with_one_pole:=proc(f,x,y,point)
	local d,j,f_infty,ib1,ib2,den1,den2,d1,d2,ansatz,a,f1,f2,equations,
	 ext,i,zero,point1,deg_inf;
	global g_conversion1,g_conversion2,`algcurves/residue`,z;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	d:=degree(f,{x,y});
	f_infty:=subs(x=1,homogeneous(f,x,y,z,polynom));
	ib1:=integral_basis(f,x,y);
	ib2:=integral_basis(f_infty,z,y,{[z,4]});
	den1:=expand(denom(ib1[d]));
	den2:=expand(denom(ib2[d]));
	if point[2]=`infinity genus1` then
		# Syntax for genus1, the function point1 may not be changed
		# so we only determine the upper bound for the denominators now
		point1:=point[1];
		# den1:=evala(Lcm(den1,denom(point1)))
		den1:=numer(g_normal(den1*denom(point1)
		 /evala(Gcm(den1,denom(point1)))))
	elif point[2]=infinity then
		deg_inf:=degree(den1,x)-degree(point[1],{x,y})+1;
		point1:=x^deg_inf*point[1]/den1;
		if deg_inf<0 then den1:=evala(Expand(den1*x^(-deg_inf))) fi
	else
		point1:=point[1]
	fi;
	# den1 is an upperbound for denominators
	d1:=degree(den1,x);
	d2:=ldegree(den2,z);
	d2:=max(d2,degree(expand(numer(point1)),{x,y})
	 -degree(expand(denom(point1)),{x,y}));
	# d2 is an upperbound for the denominators in infinity
	ansatz:=convert([seq(seq(a[i,j]*x^i*y^j,i=0..d1+d2-j)
	 ,j=0..d-1)],`+`)/den1;
	if point[2]='finite' then
		f1:=ansatz;
		f2:=ansatz+point1
	else
		f1:=ansatz+point1;
		f2:=ansatz
	fi;
	# Now we solve the following set of linear equations:
	# {f1 has no poles in 'finite', and f2 has no poles in infinity}
	ext:=g_ext([f,f1,f2]);
	f1:=g_normal(f1);
	den1:=expand(denom(f1));
	f1:=subs(g_conversion1,expand(numer(f1)));
	f2:=g_normal(subs(x=1,homogeneous(f2,x,y,z,'ratfunction')));
	den2:=expand(denom(f2));
	den2:=z^ldegree(den2,z);
	f2:=subs(g_conversion1,expand(numer(f2)));
	ib1:=subs(g_conversion1,[seq(g_expand(g_normal(i*den1),ext),i=ib1)]);
	ib2:=subs(g_conversion1,[seq(g_expand(g_normal(i*den2),ext),i=ib2)]);
	zero:=f1;
	for i from d-1 by -1 to 0 do while degree(coeff(zero,y,i),x)
	 >=degree(lcoeff(ib1[i+1],y),x) and coeff(zero,y,i)<>0 do
		zero:=g_expand(zero-ib1[i+1]*
		g_normal(lcoeff(coeff(zero,y,i),x)/lcoeff(lcoeff(ib1[i+1],y),x))
		*x^(degree(coeff(zero,y,i),x)-degree(lcoeff(ib1[i+1],y),x))
		,ext);
		if (indets([f,point],name) minus {x,y})<>{} then
			zero:=expand(collect(zero,{x,y},'distributed',g_normal))
		fi
	od od;
	equations:={a[d1,0],coeffs(zero,[x,y])};
	# These equations state that f1 is integral. We've added and 1 extra
	# linear condition, to make the solution unique. We still have to
	# add the equations for infinity:
	zero:=f2;
	for i from d-1 by -1 to 0 do while degree(coeff(zero,y,i),z)
	 >=degree(lcoeff(ib2[i+1],y),z) and coeff(zero,y,i)<>0 do
		zero:=g_expand(zero-ib2[i+1]*
		g_normal(lcoeff(coeff(zero,y,i),z)/lcoeff(lcoeff(ib2[i+1],y),z))
		*z^(degree(coeff(zero,y,i),z)-degree(lcoeff(ib2[i+1],y),z))
		,ext);
		if (indets([f,point],name) minus {x,y})<>{} then
			zero:=expand(collect(zero,{x,y},'distributed',g_normal))
		fi
	od od;
	equations:=subs(g_conversion2,{op(equations),coeffs(zero,[y,z])});
	d:=g_normal(subs(solve(equations,indets(equations,name)
	 minus indets(f)),ansatz+point1));
	if indets(d,name) minus indets([args]) <> {} then
		ERROR(`can not find parameter`,f=evala(AFactor(f)))
	fi;
	d
end:

#savelib ('iss94','function_with_one_pole','`algcurves/iss94.m`'):