Constant:=function(n,e)
#This function returns a vector all of whose n entries equal to e.
	return List([1..n], a->e);
end;

ExponentsOfGivenDegree:=function(n,d)
#This function returns all multi-indices with n entries and total 
#degree d.
	return OrderedPartitions(d+n,n)-Constant(n,1);
end;

Power:=function(a,m,one)
#This function computes the product (a_1^m_1)*(a_2^m_2)*..; both 
#lists should be equally long. The elements a_i are assumed to be 
#from a monoid, the one of which is the third parameter. The m_i 
#must be natural numbers.
local i,r;
	r:=one;
	for i in [1..Length(m)] do 
		if m[i]<>0 then r:=r*a[i]^m[i]; fi; 
	od; 
	return(r);	
end;

CoefficientsOfLinearPart:=function(u,l)
#Here, u is an element of a Universal Enveloping Algebra of a Lie
#algebra of dimension at least l. The function returns, from the
#linear part of u, only the variables with index>l, with their
#coefficients. Hence, for l=6 and 
#u=[(-1)*x.1*x.7+(1)*x.3*x.7*x.8+(1)*x.6*x.8
#	+(5)*x.6+(-1)*x.7+(2)*x.8]
#the list [[7,-1],[8,2]] is returned.
local ret, e, i;
	ret:=[];
	e:=ExtRepOfObj(u)[2];
	for i in [1..Length(e)/2] do 
		if (Length(e[2*i-1])=2) and (e[2*i-1][2]=1) and
				(e[2*i-1][1]>l)
			then AddSet(ret,[e[2*i-1][1],e[2*i]]);
		fi;
	od;
	return(ret);
end;

Blattner:=function(g,Z,n,d)
#The parameters are:
#g:	a Lie algebra
#Z:	a linear basis of g.
#n:	The elements Z[1], Z[2],...,Z[Length(Z)-n] are supposed 
#	to span a subalgebra k.
#d:	a natural number
#This function computes a realization of the pair (g,k). More
#precisely, it returns a pair (L,W), where W is the Weyl algebra 
#in n variables, and L, whose entries are elements of W, is the 
#image of Z under the realization, with coefficients truncated 
#after degree d.

local k,K,U,W,one,X,Y,x,D,L,e,exps,
Ymons,xmons,facs,i,j,coeffs,t;
	k:=Dimension(g)-n; 
		#This is the dimension of k.
	K:=LeftActingDomain(g);
		#The field.

	U:=UEA(g,Z);
		#The universal enveloping algebra of g with 
		#PBW-basis corresponding to Z.

	one:=Identity(U);
	X:=GeneratorsOfAlgebraWithOne(U);
	Y:=X{[k+1..k+n]};
		#These form a basis of g complementary to k if 
		#we identify both Lie algebras with subspaces 
		#of U.
	
	W:=WeylAlgebra(K,n);
	x:=GeneratorsOfAlgebraWithOne(W){[1..n]};
	D:=GeneratorsOfAlgebraWithOne(W){[n+1..2*n]};
		#These will be used to print the result.
	
	L:=Constant(k+n,Zero(W));

	for e in [0..d] do
		exps:=ExponentsOfGivenDegree(n,e);
		Ymons:=List(exps,m->Power(Y,m,one));
			#These are the monomials Y^m in U
		xmons:=List(exps,m->Power(x,m,Identity(W)));
			#These are the monomials x^m
		facs:=List(exps,m->Product(List(m,Factorial)));
			#This is a list of factorials m!
		for i in [1..Length(exps)] do 
			for j in [1..k+n] do
			  coeffs:=
			  CoefficientsOfLinearPart(Ymons[i]*X[j],k);
			#This computes the relevant linear part
			#of Y^m X[j]
			  for t in coeffs do 
			 	L[j]:=L[j]+
				(t[2]/facs[i])*xmons[i]*D[t[1]-k];
			  od;
			od;
		od;
	od;
	return [L,W];
end;