########################################################################################################################################
#                                                                                                               		       #
#                         We consider a first-order system of ODEs of the form                                                         #
#                                         A(x) Theta_k(y(x)) + B(x) y(x)=0,                                                            #
#      	where B(x) in K[[x]]^{n x n} and A(x) is of the form x^{alpha}Q(x) with alpha a diagonal matrix whose diagonal entries are     #	
#      		nonnegative integers in increasing ordre and Q(x) in K[[x]]^{n x n} s.t. Q(0) is the identity matrix.                  #
#        Let L= A(x) Theta_k + B(x). We present here an algorithm which computes two invertible matrices S and T                       #
#					and a k-simple operator F such that F=SLT.                                                     # 
########################################################################################################################################

			
#KSimpleForm(A,B,x,lambda,k)
#Input : A, B - The matrices of the system A*Theta_k+B with A in the set A_n. 
#	 x - the variable
#        lambda - a parameter
#	 k - a nonnegative integer
#Output : Four matrices A1, B1, S and T where S and T are two invertible matrices such that S.(A*Theta_k+B).T=A1*Theta_k+B1 is k-simple. 
#	  Det corresponds to the determinant of the matrix pencil subs(x=0,A1*lambda+B1).

KSimpleForm:=proc(A,B,x,lambda,k)

	local t0, n, A1, B1, alpha, r, i, L0, Det, S, T, q, ker_cr, ker_cr_dim, A1_0;

	t0 := time();
	n := LinearAlgebra:-RowDimension(A);

	A1 := Matrix(n,n, proc(i,j) A[i,j] end);
	B1 := Matrix(n,n, proc(i,j) B[i,j] end);

	alpha := [seq(ldegree(A1[i,i],x),i=1..n)];
	
	r := 0;   # r is the rank of A1(0)
	for i from 1 to n do
		if alpha[i]=0 then
			r := r+1;
		end if;
	end do;
	
	L0 := subs(x=0,A1)*lambda+subs(x=0,B1);
	Det := LinearAlgebra:-Determinant(L0);

	S :=  LinearAlgebra:-IdentityMatrix(n);
	T :=  LinearAlgebra:-IdentityMatrix(n);

	while Det=0 do
		A1,B1,S,T,q,ker_cr := QTCD(A1,B1,L0,x,lambda,S,T,r,n,0);
		
		T := Matrix(n,n,proc(i,j) if j<= q then x*T[i,j] else T[i,j] fi end);
		S := Matrix(n,n,proc(i,j) if i<= q then x^(-1)*S[i,j] else S[i,j] fi end);

		A1 := Matrix(n,n,proc(i,j) if i<=q and j<=q then A1[i,j] else if j>q and i<=q then simplify(x^(-1)*A1[i,j]) else if j<=q and i>q then simplify(x*A1[i,j]) else A1[i,j] fi fi fi end);
		B1 := Matrix(n,n,proc(i,j) if i<=q and j<=q then B1[i,j]+x^k*A1[i,j] else if j>q and i<=q then simplify(x^(-1)*B1[i,j]) else if j<=q and i>q then simplify(x*B1[i,j]+x^(k+1)*A1[i,j]) else B1[i,j] fi fi fi end);
		#LinearAlgebra:-Determinant(A1*lambda+B1);

		A1_0 := subs(x=0,A1);
		T := Matrix(n,n,proc(i,j) if j<= q then T[i,j] else T[i,j]-add(T[i,s]*A1_0[s,j],s=1..q) fi end);	

		A1 := Matrix(n,n,proc(i,j) if j<= q then A1[i,j] else A1[i,j]-add(A1[i,s]*A1_0[s,j],s=1..q) fi end);
		B1 := Matrix(n,n,proc(i,j) if j<= q then B1[i,j] else B1[i,j]-add(B1[i,s]*A1_0[s,j],s=1..q) fi end);	
	
		ker_cr_dim := nops(ker_cr);

		A1,B1,alpha,r,S,T := RowDependent(A1,B1,x,alpha,ker_cr,ker_cr_dim,S,T,k,r,n);		
						
		L0 := subs(x=0,A1)*lambda+subs(x=0,B1);
		Det := LinearAlgebra:-Determinant(L0);				
	end do;

	return(A1,B1,S,T,Det, time()-t0);

end proc:

# Kronecker delta
KD := proc(i,j)
	if i=j then return(1) else return(0) end if;
end proc:

#RowDependent(A,B,x,alpha,ker_cr,ker_cr_dim,S,T,k,r,n)
#Input : A, B - The matrices of the system A*Theta_k+B
#	 x - the variable
#	 alpha - The list of the diagonal entries valuations of A.
#	 ker_cr - The (left) kernel of the sub-matrix of B(0) corresponding to the constant rows of the matrix pencil
# 	 ker_cr_dim - Dimension of ker_cr
#        S, T -  The two transformations.
#	 k - a nonnegative integer 
#	 r - rank of A(0)
#	 n - size of A and B.
#Output : Matrices newA and newB defining an equivalent system to A*Theta_k+B and obtained after applying Proposition 4.3.1 ker_r_dim times. 
#	  beta contains the diagonal entries valuations of newA. 
#	  R is the rank of newA evaluated at x=0.  	  

RowDependent := proc(A,B,x,alpha,ker_cr,ker_cr_dim,S,T,k,r,n)

	local newA, newB, beta, R, ker, ker_dim, k0, i, v, i0, nu, w, p, newS, newT;
	
	ker := ker_cr;
	ker_dim := ker_cr_dim;	

	newA := Matrix(n,n, proc(i,j) A[i,j] end);
	newB := Matrix(n,n, proc(i,j) B[i,j] end);

	newS := Matrix(n,n, proc(i,j) S[i,j] end);
	newT := Matrix(n,n, proc(i,j) T[i,j] end);
	
	beta := alpha;
	R := r;

	ker := convert([op(ker)], Matrix);
	ker :=  LinearAlgebra:-GaussianElimination( LinearAlgebra:-Transpose(ker));

	k0:=1;

	for i from 1 to ker_dim do
		
		v := LinearAlgebra:-Row(ker,i);
		while v[k0]=0 do
			k0 := k0+1;
		end do;
		v := map(simplify,1/v[k0]*v);
		i0 := r+k0;
		
		newS := Matrix(n,n, proc(i,j) if i<>i0 then newS[i,j] else add(v[k0+s]*newS[i0+s,j],s=0..n-i0) fi end);
		newA := Matrix(n,n, proc(i,j) if i<>i0 then newA[i,j] else add(v[k0+s]*newA[i0+s,j],s=0..n-i0) fi end);
		newB := Matrix(n,n, proc(i,j) if i<>i0 then newB[i,j] else add(v[k0+s]*newB[i0+s,j],s=0..n-i0) fi end);		

		newT := Matrix(n,n, proc(i,j) if j<=i0 then newT[i,j] else simplify(-KD(beta[i0],beta[j])*v[j-r]*newT[i,i0]+newT[i,j]) fi end);
		newA := Matrix(n,n, proc(i,j) if j<=i0 then newA[i,j] else simplify(-KD(beta[i0],beta[j])*v[j-r]*newA[i,i0]+newA[i,j]) fi end);
		newB := Matrix(n,n, proc(i,j) if j<=i0 then newB[i,j] else simplify(-KD(beta[i0],beta[j])*v[j-r]*newB[i,i0]+newB[i,j]) fi end);

		nu := min(map(ldegree, LinearAlgebra:-Row(newB,i0),x), degree(newA[i0,i0],x));

		beta[i0] := beta[i0]-nu;

		if nu = degree(newA[i0,i0],x) then
			R:=R+1;
		end if;

		newS := Matrix(n,n, proc(i,j) if i=i0 then simplify(x^(-nu)*newS[i,j]) else newS[i,j] fi end);
		newA := Matrix(n,n, proc(i,j) if i=i0 then simplify(x^(-nu)*newA[i,j]) else newA[i,j] fi end);
		newB := Matrix(n,n, proc(i,j) if i=i0 then simplify(x^(-nu)*newB[i,j]) else newB[i,j] fi end);
	end do;

	w := [seq(i=beta[i],i=1..n)];
  	p := proc(a,b) rhs(a) < rhs(b) end:
  	w := sort(w,p);
	
	beta := [seq(rhs(w[i]),i=1..n)];
	newB := Matrix(n,n, proc(i,j) newB[lhs(w[i]),j] end);
	newB := Matrix(n,n, proc(i,j) newB[i,lhs(w[j])] end);

	newA :=Matrix(n,n, proc(i,j) newA[lhs(w[i]),j] end);
	newA := Matrix(n,n, proc(i,j) newA[i,lhs(w[j])] end);

	newS := Matrix(n,n, proc(i,j) newS[lhs(w[i]),j] end);
	newT := Matrix(n,n, proc(i,j) newT[i,lhs(w[j])] end);
	
	return(newA,newB,beta,R,newS,newT);

end proc:

#QTCD(A,B,MP,x,lambda,S,T,r,n,q)
#Input : A, B - The matrices of the system A*Theta_k+B
#        MP - the matrix pencil associated with the system
#	 x - the variable of the system
#        lambda - a parameter
#	 alpha - the list of the diagonal entries valuations of A 
#	 S,T - The two transformations
#	 r - the rank of A(0)
#	 n - the size of the system
#	 q - an integer
#Output: Two matrices newA and newB such that the matrix pencil of the system newA*Theta_k+newB 
#	 is in the form (4.15) and satisfy condition (4.16).

QTCD:=proc(A,B,MP,x,lambda,S,T,r,n,q)

	local newA, newB, newS, newT, Sub_L0, new_q, CR, ker_cr, ker, v, k0, a, L00;

	new_q := q;

	newA := Matrix(n,n, proc(i,j) A[i,j] end);
	newB := Matrix(n,n, proc(i,j) B[i,j] end);

	newS := Matrix(n,n, proc(i,j) S[i,j] end);
	newT := Matrix(n,n, proc(i,j) T[i,j] end);

	CR := LinearAlgebra:-SubMatrix(MP,[r+1..n],[new_q+1..n]);
	ker_cr := LinearAlgebra:-NullSpace(LinearAlgebra:-Transpose(CR));
	
	if nops(ker_cr) <> 0 then 
		return(newA,newB,newS,newT,new_q,ker_cr);
	else
		Sub_L0 := LinearAlgebra:-SubMatrix(Matrix(n,n, proc(i,j) MP[i,j] end),[new_q+1..n],[new_q+1..n]);
		L00 := subs(lambda=0,Sub_L0);
		ker := LinearAlgebra:-NullSpace(LinearAlgebra:-Transpose(L00));
		
		v := ker[1];
	
		k0 := 1;
	
		while v[k0]=0 do
			k0 := k0+1;
		end do;
		
		if k0<>1 then
			a := v[k0];
			v[k0] := v[1];
			v[1] := a;
		end if;
	
		v := 1/v[1]*v;

		newS := LinearAlgebra:-RowOperation(newS,[new_q+1,k0+new_q]); 
		newA := LinearAlgebra:-RowOperation(newA,[new_q+1,k0+new_q]); 
		newB := LinearAlgebra:-RowOperation(newB,[new_q+1,k0+new_q]); 
		
		newS := Matrix(n,n, proc(i,j) if i<>new_q+1 then newS[i,j] else add(v[s-new_q]*newS[s,j],s=new_q+1..n) fi end);
		newA := Matrix(n,n, proc(i,j) if i<>new_q+1 then newA[i,j] else add(v[s-new_q]*newA[s,j],s=new_q+1..n) fi end);
		newB := Matrix(n,n, proc(i,j) if i<>new_q+1 then newB[i,j] else add(v[s-new_q]*newB[s,j],s=new_q+1..n) fi end);
		
		newT := LinearAlgebra:-ColumnOperation(newT,[new_q+1,k0+new_q]);
		newA := LinearAlgebra:-ColumnOperation(newA,[new_q+1,k0+new_q]);	
		newB := LinearAlgebra:-ColumnOperation(newB,[new_q+1,k0+new_q]);		
		
		newT := Matrix(n,n, proc(i,j) if j <= new_q+1 or j>r then newT[i,j] else simplify(-v[j-new_q]*newT[i,new_q+1]+newT[i,j]) fi end); 
		newA := Matrix(n,n, proc(i,j) if j <= new_q+1 or j>r then newA[i,j] else simplify(-v[j-new_q]*newA[i,new_q+1]+newA[i,j]) fi end); 
		newB := Matrix(n,n, proc(i,j) if j <= new_q+1 or j>r then newB[i,j] else simplify(-v[j-new_q]*newB[i,new_q+1]+newB[i,j]) fi end); 

		new_q := new_q+1;
		 
		QTCD(newA,newB,subs(x=0,newA*lambda+newB),x,lambda,newS,newT,r,n,new_q);
	end if;
end proc: