with (powseries):
with(LinearAlgebra):
with (linalg):
with(PolynomialTools);
unprotect(normalize);

ExpandsF := proc(VV,N) local e;
e:= evalpow(VV);
convert(tpsform(e,z,N),polynom)
end:
"ExpandsF(function_of_z, to order N)";


 CanPolysV := proc(VV,N) local w,n,g,i,W,m; global U,v;
	unassign('U','v');
	W:=simplify(1/diff(VV,z));
 	n:=0;
	w[0]:=coeftayl(W,z=0,0);
	g[0]:=1;
                        
	 for i from 1 to N do:
 		n:=n+1;w[n]:=coeftayl(W,z=0,n);  
 		g[n]:=expand(x*simplify(w[0]*g[n-1]+sum(w[j]*diff(g[n-1],x$j),j=1..n))):print(g[n]);
 	end do:
		print(`U=`);
		U := add(v^m/factorial(m)*coeff(g[m],x),m=0..n)
end:
	"CanPolysV(function_of_z_to_invert, number of polys to generate)";

 CanPolysW := proc(VV,N) local w,n,g,i,W,m; global U,v;
	unassign('U','v');
	W:=VV;
 	n:=0;
	w[0]:=coeftayl(W,z=0,0);
	g[0]:=1;
                        
 	for i from 1 to N do:
 		n:=n+1;w[n]:=coeftayl(W,z=0,n);  
 		g[n]:=expand(x*simplify(w[0]*g[n-1]+sum(w[j]*diff(g[n-1],x$j),j=1..n))):print(g[n]);
 	end do:
		print("U=");
		 U := add(v^m/factorial(m)*coeff(g[m],x),m=0..n)
end:"CanPolysW(W = 1 over V', number of polys to generate)";

kron:= proc(AA,BB) local i,j,CX,UX,ii,jj,n,m;
n:=rowdim(AA);m:=coldim(AA);
for i to n do: for j to m do:
CX[j]:=scalarmul(BB,AA[i,j])
od;
UX[i]:=concat(seq(CX[jj],jj=1..m))
od;
stackmatrix(seq(UX[ii],ii=1..n))
end:"kron(A,B) returns the Kronecker product of two matrices";

cmm:=proc(XX,YY)
    evalm(XX&*YY-YY&*XX)
end:"cmm(X,Y) returns the commutator of two matrices";

normalize:=proc(f) local s,kk,t,i,d,zz,of,cf;
s:=0:
cf:=[coeffs(f,[seq(z||k,k=1..N)],'of')];

	for i to nops(cf) do 
		for kk to N do 

		d[kk]:=degree([of][i],z||kk);	

	if (d[kk]=0) 
			then zz[kk]:=JJ 
			else zz[kk]:=(z||kk)^(d[kk]) 
		fi 
			od; 

 	       s:=s+cf[i]*`&*`(seq(zz[a],a=1..N)) 
			       od;
	s
end:"normalize converts multiplication into matrix multiplication";

invbymatW:=proc(ff,order) local f,dg,CL,W,M,Q,MW,UMAT,VMAT,vv,xx,EXU,i,k;
	unassign('z','v','x'):
	f:=convert(taylor(ff,z=0,order+1),polynom);
	dg:=degree(f,z);

	CL:=CoefficientList(f,z);

	if(dg<order) then
		for k to (order-dg) do 
			CL:=[op(CL),0] 
		od;
	fi;		

	dg:=order;

	W:=matrix(dg,dg,(i,j)->if(i-j>=-1) then CL[i-j+2] else 0 fi);
	M:=diag(seq(k,k=1..dg));	
	Q:=diag(seq(1/GAMMA(k),k=1..dg));
	MW:=multiply(M,W);

		 UMAT[1]:=Vector(dg,i->if(i=1) then CL[1] else 0 fi);

	 for i from 2 to dg do 
		UMAT[i]:=multiply(UMAT[i-1],MW);
	od:	 	
		vv:=vector([seq(v^i/i!,i=1..dg)]);
		xx:=vector([seq(x^k,k=1..dg)]);
	 	EXU:=multiply(vv,multiply(stackmatrix(seq(UMAT[k],k=1..dg)),Q,xx));

		EXU, coeff(EXU,x)
end:"invbymatW(W,order) outputs exponential of xU(v) and U(v) to that order in v";
	
invbymatV:=proc(ff,order) local WW;
	unassign('z'):
	invbymatW(1/diff(ff,z),order)
end:"invbymatV(V,order) outputs exponential of xU(v) and U(v) to that order in v";