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

macro(
	genus1=`algcurves/genus1`,
	genus1_compute_x=`algcurves/g1_comp_x`
):

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`
):

# Input: a polynomial f in x and y with genus 1, the variables x,y,x0,y0
#  and a regular point on the curve.
#  f must be irreducible over Qbar
# If the last argument is `no inverse` then the inverse of the isomorphism
#  is not computed (i.e. the image of x and y is not computed).
# If [0,1,0] is a point on the curve then:
#   - with the option `no inverse` the output is not necessarily
#     a polynomial in y.
#   - without the option `no inverse` it does become a polynomial
#     by doing a division in the function field
# Output: a list containing:
#  f0 such that Qbar(x)[y]/(f) is isomorphic with Qbar(x0)[y0]/(f0)
#  The image of x0 under this isomorphism
#  The image of y0
#  The image of x under the inverse isomorphism
#  The image of y
# f0 is of the form y0^2 + polynomial of degree 3 in x0^2
genus1:=proc(ff,x,y,x0,y0,point)
	global `algcurves/residue`;
	local a,d,t,v,i,f,f0,x0v,y0v,ansatz,p2,p3,p;
	options remember, 
	    `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	f:=collect(ff,{x,y},'distributed',g_normal);
#if nargs<6 or not type(point,list) then
#	# No point was specified, compute one:
#	RETURN(procname(args[1..5],find_a_point(f,x,y),args[6..nargs]))
if not member(`no inverse`,{args}) then
	v:=procname(
	 # remove the `j invariant` argument, otherwise the options remember
	 # in this procedure does not help.
	 `if`(args[nargs]=`j invariant`,args[1..nargs-1],args)
	 ,`no inverse`);
	f0:=v[1];
	if member(`j invariant`,{args}) then
		f0:=subs(y0=0,f0);
		f0:=f0/lcoeff(f0,x0);
		a:=coeff(f0,x0,1);
		d:=coeff(f0,x0,0);
		RETURN(evala(6912*a^3/(27*d^2+4*a^3)))
	fi;
	# Remove possible y's in the denominators:
	v:=seq(FFDIV(f,x,y,v[i]),i=2..3);
	# v=x0v,y0v
	RETURN([f0,v,genus1_compute_x(v,x,y,f,f0,x0,y0),
	 genus1_compute_x(v,y,x,f,f0,x0,y0)])
fi;
	d:=degree(f,{x,y});
	if coeff(f,y,d)=0 then
	  if coeff(f,x,d)=0 then
		# [0,1,0] is a point on the curve, we want to avoid this
		v:=[point[1]-point[2],point[2],point[3]],args[7..nargs];
		# apply transformation on f to remove the point [0,1,0]
		v:=procname(subs(x=x+y,f),args[2..5],v);
		if nops(v)>1 then
			# Do the reverse transformation
			i:=subs(x=x-y,[v[2],v[3]]);
			if nops(v)>3 then
				v:=[v[1],op(i),v[4]+v[5],v[5]]
			else
				v:=[v[1],op(i)]
			fi
		fi;
		RETURN(v)
	  else
		# [0,1,0] is a point on the curve, we want to avoid this
		v:=[point[2],point[1],point[3]],args[7..nargs];
		v:=procname(subs({x=y,y=x},f),args[2..5],v);
		if nops(v)>1 then
			# Do the reverse transformation
			i:=subs({x=y,y=x},[v[2],v[3]]);
			if nops(v)>3 then
				v:=[v[1],op(i),v[5],v[4]]
			else
				v:=[v[1],op(i)]
			fi
		fi;
		RETURN(v)
	  fi
	fi;
	# Now coeff(f,y,d)<>0, so [0,1,0] is not a point on the curve
	# Compute a function with a pole in the given point:
	if point[3]=0 then # point in infinity
		p:=[(homogeneous(evala(Quo(subs({t=0,x=1}
		 ,homogeneous(f,x,y,t,polynom)),y-point[2]/point[1],y))
		 ,y,t,x,polynom))/x^(degree(f,y)-2),`infinity genus1`]
	else
		p:=[evala(Quo(evala(Expand(subs(x=point[1]/point[3],f)))
		 ,y-point[2]/point[3],y))/(x-point[1]/point[3]),'finite']
	fi;
	# Now make it a double pole with an indeterminate as residue
	p2:=[evala(Normal(rem(expand(-p[1]^2+`algcurves/residue`
	 *p[1]),f,y))),p[2]];
	# Compute a function with this pole and no other poles:
	x0v:=function_with_one_pole(f,x,y,p2);
	p3:=[evala(Normal(rem(expand(-p[1]*x0v+`algcurves/residue`
	*p[1]),f,y))),p[2]];
	# Now compute a function with pole order 3:
	y0v:=function_with_one_pole(f,x,y,p3);
	# look for a relation f0 between x0v and y0v using an ansatz:
	i:=0;
	ansatz:=y0^2+x0^3+a[1]*y0+a[2]*x0*y0+a[3]+a[4]*x0+a[5]*x0^2;
	while has(ansatz,{a[1],a[2],a[3],a[4],a[5]}) do
		i:=i+1;
		# avoid certain values i for x
		while g_normal(subs(x=i,denom(x0v)*denom(y0v)))=0 do
			i:=i+1
		od;
		# find linear equations in a[1] .. a[5]
		ansatz:=subs(solve({coeffs(expand(rem(expand(numer(
		 subs(x0=x0v,y0=y0v,x=i,ansatz))),expand(subs(x=i,f)),y)),y)}
		 ,{a[1],a[2],a[3],a[4],a[5]}),ansatz)
	od; # repeat until all a[i] are determined (usually in the 1'st step)
	f0:=ansatz;
	# Now normalize f0 a little further:
	y0v:=g_normal(y0v+subs(x0=x0v,coeff(f0,y0,1))/2);
	# clear the coefficient of y0^1:
	f0:=collect(subs(y0=y0-coeff(f0,y0,1)/2,f0),
	 [x0,y0],'distributed',g_normal);
	x0v:=g_normal(x0v+coeff(f0,x0,2)/3);
	# clear the coefficient of x0^2:
	f0:=collect(subs(x0=x0-coeff(f0,x0,2)/3,f0),
	 [x0,y0],'distributed',g_normal);
	if has('Weierstrass',{args}) then # Weierstrass normal form
		x0v:=-x0v;
		y0v:=2*y0v;
		f0:=y0^2+4*subs(x0=-x0,y0=0,f0)
	fi;
	[f0,x0v,y0v]
end:

# Express x as an expression in x0 and y0.
# f0 = algebraic relation between x0 and y0
# f = algebraic relation between x and y
# x0v = image of x0 in Qbar(x,y)
# y0v = image of x0 in Qbar(x,y)
genus1_compute_x:=proc(x0v,y0v,x,y,f,f0,x0,y0)
	local cp,result,v,i;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	cp:=genus1_charpol(x0v,x,y,f,x0);
	result:=-coeff(cp,x,1)/2;
	cp:=g_normal((x^2-subs(x=x+result,cp))/coeff(f0,y0,0));
if cp<>0 then
	# Now add y0*sqrt(cp) to result:
	v:=[seq(evala(Sqrfree(i),'expanded')[2],i=[numer(cp),denom(cp)])];
	v:=g_normal(convert([seq(i[1]^(i[2]/2),i=v[1])],`*`)/
		    convert([seq(i[1]^(i[2]/2),i=v[2])],`*`));
	cp:=factors(numer(x^2+g_normal(cp/v^2))
	 ,indets([args],RootOf))[2];
	if has(cp,x0) then bug() fi;
	for i in cp do if has(i,x) and degree(i[1],x)=1 then
		v:=v*subs(solve({i[1]},{x}),x)*y0;
		result:=result+v,result-v;
		break
	fi od;
fi;
	i:=0;
	result:={result};
	while nops(result)>1 do for v in result do
		i:=i+1;
		# avoid certain values i for x
		while g_normal(subs(x=i,denom(x0v)*denom(y0v)))=0 do
			i:=i+1
		od;
		if g_normal(rem(expand(numer(subs(x0=x0v,y0=y0v,x=i,numer(v))
		 -i*subs(x0=x0v,y0=y0v,x=i,denom(v)))),subs(x=i,f),y))<>0 then
			result:=result minus {v};
			break
		fi
	od od;
	op(result)
end:

#savelib ('genus1','genus1_compute_x','`algcurves/genus1.m`'):