macro(
	IntBasis = algcurves['integral_basis'],
	ScalarMult = LinearAlgebra:-MatrixScalarMultiply,
	SolveLinear = SolveTools:-Linear,
	g_conversion1 = `algcurves/g_conversion1`,
	g_conversion2 = `algcurves/g_conversion2`,
	g_normal=`algcurves/g_normal`,
	g_ext = `algcurves/g_ext`,
	truncate_subs = `algcurves/truncate_subs`,
	ext_to_coeffs=`algcurves/e_to_coeff`,
	lift_exp_m1 = `algcurves/lift_exp_m1`,
	`puiseux/technical_answer`=`algcurves/puiseux_te`,
	RAD={'radical','nonreal'}
):


NumerDenom := proc(L,x,y)
	local d,B,v,i;
	# options trace;
	d := 1;
	B := L;
	v := select(has,map(denom,B),x);
	while v<>[] do
		v := v[-1];
		B := [seq(Normalizer(i*v),i=B)];
		d := d*v;
		v := select(has,map(denom,B),x);
	od;
	map(collect,B,y), d
end:

# Express Bn as linear combinations of Bi.
ExpressAsLinComb := proc(Bn, Bi, x, y, n)
	local B,d,M,i,pol,j,c,cd,k;
	# options trace;
	B := Bn;
	d := 1;
	M := Matrix(n,n,0);
	for i to n do
		pol := B[i];
		for j from i to 1 by -1 do
			c := Normalizer(coeff(pol,y,j-1)/lcoeff(Bi[j],y));
			cd := denom(c);
			if has(cd,x) then
				d := cd*d;
				B := [seq(collect(k*cd,y),k=B)];
				M := ScalarMult(M,cd);
				pol := collect(cd*pol,y);
				c := numer(c)
			fi;
			M[i,j] := c;
			pol := collect(pol-c*Bi[j],y)
		od
	od;
	map(collect,M,x,Normalizer), d
end:


# v is a triangular basis of some C[x]-module
#    (triangular means that degree(v[i],y) = i-1).
# a is a C(x)-linear combination of v.
# Now compute a triangular basis of the C[x]-module C[x]*v + C[x]*a.
InsertInBasis := proc(a, v, x, y)
	local d,l,q,A,R,s,t;
	if a=0 then return v fi;
	d := degree(a,y);
	l := Normalizer(lcoeff(a,y));
	if l=0 then return procname(collect(a,y,Normalizer),v,x,y) fi;
	q := Normalizer(l/lcoeff(v[d+1],y));
	gcdex(numer(q), denom(q), 1, x, 's', 't');
	A := collect(s*a+t*v[d+1],y,Normalizer);
	R := Normalizer(l/lcoeff(A,y));
	if has(denom(R),x) then
		error "division failed"
	fi;
	R := collect(a - R*A,y,Normalizer);
	procname(R, [op(1..d,v),A,op(d+2..-1,v)], x, y)
end:


degreesList := proc(M, n, x)
	local i,j;
	[seq(max(seq(degree(M[i,j],x),j=1..n)),i=1..n)]
end:

TransposeLcMatrix := proc(M, n, x, D)
	Matrix(n,n,[seq([seq(coeff(M[i,j],x,D[i]),i=1..n)],j=1..n)])
end:

NormalBasis := proc(M, n, x)
	local P,PM,D,N,v,s,i,j,vP,vPM;
	P := LinearAlgebra:-IdentityMatrix(n,'compact'=false);
	PM := Matrix(M);
	do
		D := degreesList(PM, n, x);
		N := TransposeLcMatrix(PM, n, x, D);
		v := LinearAlgebra:-NullSpace(N);
		if v={} then break fi;
		v := v[1];
		for s from n by -1 do if v[s]<>0 then break fi od;
		for j from n-1 to 1 by -1 do if D[j]>D[s] and v[j]<>0 then
			s := j
		fi od;
		v := Vector['row']([seq(Normalizer(v[j]/v[s])*x^(D[s]-D[j]),j=1..n)]);
		vP := v . P;
		vPM := v . PM;
		for j to n do
			P[s,j] := collect(vP[j],x); # ,Normalizer);
			PM[s,j] := collect(vPM[j],x); # ,Normalizer)
		od;
	od;
	P, [seq(D[i]-D[1],i=1..n)]
end:



PolarPartsAtPlaces := proc(F, f, x, y, ext)
	local v,i,ex,f1,F1,dF,nF,res,p,r,a,P,E;
	if indets([args],RAD)<>{} then
		v := [args];
		i := radfield(indets(v,RAD));
		eval(procname(op(eval(v,i[1]))),i[2])
	elif coeff(lcoeff(f,y),x,0) = 0 then
		procname(subs(y=y/x,F), primpart(subs(y=y/x,f),y), x,y,ext)
	else
	ex := g_ext([args]);
	f1 := subs(g_conversion1,f);
	v := `puiseux/technical_answer`(f1,x,y,0,ex);
	F1 := Normalizer(F);
	dF := ldegree(collect(denom(F1),x,Normalizer),x);
	nF := collect(numer(F1),[y,x],Normalizer);
	res := NULL;
	for p in v do
		r := degree(p[3],x);
		# Determine the required accuracy:
		a := dF - min(seq(ldegree(coeff(nF,y,i),x)
		 +(i-1)*ldegree(p[1]+co[1]*x^p[2],x)/r,i = 1..degree(nF,y)));
		# Lift to this accuracy:
		P := lift_exp_m1(p,f,x,y,a);
		E := truncate_subs(subs(g_conversion1, x=P[3], nF),
		 x, y, P[1], r*dF, P[4]);
		res := res, [seq( # `if`(coeff(E,x,i)=0,NULL,
		 {ext_to_coeffs(coeff(E,x,i),ex)},i=ldegree(E,x)..r*dF-1)]
	od;
	subs(g_conversion2,{res});
	fi
end:


Candidates := proc(f,x,y)
	local n,finf,Binf,B,M,d,P,Degs,v,i,j,N,vars,ext,V,so,W,s;
	n := degree(f,y);
	finf := primpart(subs(x=1/x,f),y);
	# Integral basis at x=infinity
	Binf := IntBasis(finf,x,y,{x});
		# if location=infinity then
		#	Binf := InsertInBasis(collect(subs(x=1/x,f1),y,Normalizer),
		#	  Binf, x, y)
		# fi;
	Binf := map(collect, subs(x=1/x,Binf), y, Normalizer);

	# Integral basis for finite x
	B := IntBasis(f,x,y);
		# if location="finite" then
		#	B := InsertInBasis(f1, B, x, y)
		# fi;
		#
		# Bn, Bd := NumerDenom(B,x,y);
	M, d := ExpressAsLinComb(B, Binf, x, y,n):
	P, Degs := NormalBasis(M, n, x):
	v := [x, seq(`if`(Degs[i]=1,collect(add(P[i,j]*B[j],j=1..n),y,
	 Normalizer),NULL),i=2..n)]:
	N := nops(v);
	if N=1 then
		# Didn't find a good simplification.
		return FAIL
	fi;
	vars := {seq(co[i],i=1..N)};
	v := add(co[i]*v[i],i=1..N);
	# Format: For each Place there is a list of coefficients of the
	# expansion at that place: V = {[c11,c12,..], [c21,c22,..], ... }
	# Each c.ij itself is also a list of coeffs w.r.t. the algebraic
	# extension of that Place.
	ext := [];
	V := PolarPartsAtPlaces(subs(x=1/x,v), finf, x, y, ext);
	so := {};
	while V<>{} do
		V := subs(so,V) minus {[]};
		W := sort([op(V)],proc(a,b) evalb(nops(a[1])>nops(b[1])) end);
		for i in W do
			s := SolveLinear(subs(so,i[1]),indets(v,'name')
			 intersect vars);
			V := V minus {i};
			if map(rhs,s)<>{0} then
				v := subs(s,v);
				V := subs(s, V union {i[2..-1]});
				so := so,s
			fi
		od
	od;
	v := eval(v);
	V := [seq(coeff(v,i), i = indets(v,'name') intersect vars)];
end: