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

macro(
	genus1=`algcurves/genus1`,
	genus1_compute_x=`algcurves/g1_comp_x`,
	solve=`solve/linear`,
	function_with_one_pole=`algcurves/f_with_1_p`,
	homogeneous=`algcurves/homogeneous`,
	genus1_charpol=`algcurves/g1_charpol`,
	FFDIV=`algcurves/FFDIV`
):

# 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(f,x::name,y::name,x0::name,y0::name,point::list)
	global `algcurves/residue`;
	local a,d,t,v,i,f0,x0v,y0v,ansatz,p2,p3,p;
	options remember, 
`Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;

if not member(`no inverse`,{args}) then
		# The option `j invariant` is obsolete.
		# v:=procname(
		# `if`(args[nargs]=`j invariant`,args[1..nargs-1],args)
		# ,`no inverse`);
	v:=procname(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 Normalizer(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
	[f0,v,genus1_compute_x(v,x,y,f,f0,x0,y0),
	 genus1_compute_x(v,y,x,f,f0,x0,y0)]
else
	d:=degree(f,{x,y});
    if d>3 and member(2,{degree(f,x),degree(f,y),d-ldegree(f,{x,y})})
    and not member(`no g2`,[args]) then
	v:=`algcurves/genus2`(f,x,y,a,i,1,[],`no inverse`);
	t:=`algcurves/Weierstrassform`(v[1],a,i,x,y,args[7..nargs],`no g2`);
	[subs(x=x0,y=y0,t[1])
	 ,op(map(Normalizer,subs(a=v[2],i=v[3],[t[2],t[3]])))]
    else
	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:=[Normalizer(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:=[Normalizer(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 Normalizer(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:=Normalizer(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',Normalizer);
	x0v:=Normalizer(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',Normalizer);
	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]
    fi
fi
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,z;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`;
	cp:=genus1_charpol(x0v,x,y,f,x0);
	result:=-coeff(cp,x,1)/2;
	cp:=Normalizer((x^2-subs(x=x+result,cp))/coeff(f0,y0,0));
if cp=0 then
	result
else
	# Now add y0*sqrt(cp) to result:
	v:=[seq(sqrfree(i,x0)[2],i=[numer(cp),denom(cp)])];
	v:=Normalizer(mul(i[1]^iquo(i[2],2),i=v[1])/
		      mul(i[1]^iquo(i[2],2),i=v[2]));
	cp:=factors(numer(x^2+Normalizer(cp/v^2))
	 ,indets([args],RootOf))[2];
	if has(cp,x0) then error "should not happen" 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;
	i:=lcm(denom(x0v),denom(y0v));
	cp:=x0=Normalizer(x0v*i), y0=Normalizer(y0v*i), z=i;
	i:=0;
	result:=map(homogeneous,{result},x0,y0,z,'ratfunction');
	while nops(result)>1 do for v in result do
		userinfo(5,'algcurves',`zero testing...`);
		if Normalizer(rem(evala(Expand(
		 subs(cp,x=i,numer(v)-i*denom(v)) )),subs(x=i,f),y))<>0 then
			result:=result minus {v};
			break
		fi;
		i:=i+1
	od od;
	subs(z=1,op(result))
fi
end:

#savelib ('genus1','genus1_compute_x'):