# This file implements "Normalization of an integral basis at infinity" from:
#
# 	Chapter 2, Section 3 in:  B. Trager, Integration of algebraic functions,
#	Ph.D. thesis, Dept. of EECS, MIT, (1984)  dspace.mit.edu/handle/1721.1/15391
#
# Normalization is used in Trager's thesis and in later papers to compute
# Riemann-Roch spaces. For example, if the D is the polar divisor of x^d then L(D),
# the Riemann-Roch space for D, is computed with the program UpToDegree below.
#
# Mark van Hoeij


# Compute a basis of all h in C(x)[y]/(f) that are integral over C[x]
# for which h/x^d is integral at x=infinity.
#
UpToDegree := proc(f, x, y, d)
	local v, Degs, i, j;
	v, Degs := NormalBasis(f,x,y);
	[seq(seq(v[i]*x^j, j=0..d-Degs[i]), i=1..nops(v))]
end:


NormalBasis := proc(f,x,y) local n,finf,Binf,B,M,P,Degs,i,j;
	n := degree(f,y);
	
	# Compute Binf := Integral Basis at x=infinity
	finf := primpart(subs(x=1/x,f),y);
	Binf := algcurves['integral_basis'](finf,x,y,{x});
	Binf := map(collect, subs(x=1/x,Binf), y, Normalizer);

	# Write Integral basis for finite x in terms of Binf:
	B := algcurves['integral_basis'](f,x,y);
	M := ExpressAsLinComb(B, Binf, x, y,n)[1];

	P, Degs := NormalBasisDoIt(M, n, x):
	[seq(normal(add(P[i,j]*B[j],j=1..n)),i=1..n)], Degs
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;
	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 := LinearAlgebra:-MatrixScalarMultiply(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:

NormalBasisDoIt := 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);
			PM[s,j] := collect(vPM[j],x);
		od
	od;
	P, [seq(D[i]-D[1],i=1..n)]
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) local i,j;
	Matrix(n,n,[seq([seq(coeff(M[i,j],x,D[i]),i=1..n)],j=1..n)])
end: