#########################################################################################################################################################
#            																		#
#    		 This algorithm computes a basis of the regular formal solution space up to a fixed order (given in input)	 	   		#
#				of a simple higher-order linear differential system of the form 							#
#			   A_l(x) \vartheta^l(y(x)) + ... + A_1(x) \vartheta(y(x)) + A_0(x) y(x) =0							#
#    		 where \vartheta is the Euler derivation (\vartheta=x.d/dx) and the A_i(x)'s are square matrices of analytic functions at x=0. 		#
#																			#				
#########################################################################################################################################################     

with(LinearAlgebra):

#BCE(Liste_of_Matrices,x,set_parameters,ord)
#INPUT:		Liste_of_Matrices- The list of the coefficient matrices [A_l(x),...,A_0(x)] of the differential system. 
# 		x- variable
#		set_parameters- The set of parameters containing in the system
#		ord- the order up to which we want to compute the regular solutions 
#OUTPUT: 	The general regular solution up to order ord.

BCE:=proc(Liste_of_Matrices,x,set_parameters,ord)
	global c;
	local list, i, t0, s, y, n, l, L, InvL, p, Det_L, Deg_Det_L, succdiff, fact, count, K, Ens, d, j, dif, MPS, succdiffMPS, card_Ens_s, U_0, list_int_diff, alg_mult, list_alg_mult, max_kappa, eigenval, ind_U_0, sol, Y, nb_par; 
	t0:=time();
	c:='c';
	nb_par:=0;
	y:=[];
	n:=RowDimension(Liste_of_Matrices[1]);
	l:=nops(Liste_of_Matrices)-1;
	list:=[seq(map(series,Liste_of_Matrices[i],x=0,ord+1),i=1..l+1)]; 
	L:=add(ScalarMultiply(map(coeff,list[i],x,0), (_mu)^(l-i+1)), i = 1 .. l+1);
	InvL:=MatrixInverse(L);
	p:=max(map(degree,L,_mu));
	Det_L:=Determinant(L);
	fact:=factors(Det_L);
	Deg_Det_L:=degree(Det_L,_mu);
	count:=1;
	K:={};
	if nops(set_parameters)<>0 then
		K:={op(select(has,fact[2],_mu))};
	else
		K:={op(fact[2])};
	end if;
	while nops(K)<> 0 do
		Ens[count]:={[0,K[1]]}; 
		d:=degree(K[1][1],_mu);
		for j from 2 to nops(K) do		
			if degree(K[j][1])=d then
				#dif:=(-coeff(K[j][1],_mu,d-1)/coeff(K[j][1],_mu,d)+coeff(K[1][1],_mu,d-1)/coeff(K[1][1],_mu,d))/d;
				dif:=(coeff(K[1][1],_mu,d-1)-coeff(K[j][1],_mu,d-1))/d;
				#if type(dif,integer)=true and dif<>0 then
				if type(dif,integer)=true and K[1][1]/expand(subs(_mu=_mu+dif,K[j][1]))=1 then
					Ens[count]:={op(Ens[count]),[dif,K[j]]};
				end if;
			end if;
		end do;
		K:=K minus {seq(Ens[count][k][2],k=1..nops(Ens[count]))};
		Ens[count]:={seq([Ens[count][k][1],[RootOf(Ens[count][k][2][1],_mu),Ens[count][k][2][2]]],k=1..nops(Ens[count]))};
		count:=count+1;
	end do;  
	succdiff:=[L,seq(simplify(ScalarMultiply(map(diff,L,_mu$k),1/(k!))), k=1..min(p,Deg_Det_L-1)),seq(ZeroMatrix(n),i=1..Deg_Det_L-min(p,Deg_Det_L-1)-1)];                                                                     
	MPS:=MPSupp(list,x,n,l,ord);                                                                            
	succdiffMPS:=[seq([MPS[i],seq(ScalarMultiply(map(diff,MPS[i], _mu$k),1/(k!)), k = 1 .. Deg_Det_L-1)], i=1..ord)];                                                                                       
	for s from 1 to count-1 do
		card_Ens_s:=nops(Ens[s]);
		eigenval:=Ens[s][1][2][1];
		alg_mult:=Ens[s][1][2][2];
		sol:=LinearSolve_with_LU([seq(convert(subs(_mu=eigenval,succdiff[k]),Matrix),k=1..alg_mult)],[seq(Vector(n),i=1..alg_mult)],n,set_parameters,nb_par);
		U_0:=sol[1];
		while Equal(U_0[1],Vector(n))=true do
			U_0:=[seq(U_0[i],i=2..nops(U_0))];
		end do;
		ind_U_0:=indets(U_0) minus set_parameters;
		nb_par:=sol[2];
		max_kappa:=nops(U_0);      # the maximum value of partial multiplicity associated to eigenval
		if card_Ens_s=1 then 
			y:=[op(y),Non_Integer_Difference(x,U_0,ord,max_kappa,eigenval,succdiff,succdiffMPS,InvL)];  
		else
			list_int_diff:=[seq(Ens[s][j][1]-Ens[s][1][1],j=2..card_Ens_s)];
			list_alg_mult:=[seq(Ens[s][j][2][2],j=2..card_Ens_s)];
			Y:=Integer_Difference(x,U_0,ord,max_kappa,n,p,eigenval,succdiff,succdiffMPS,list_int_diff,list_alg_mult,InvL,set_parameters,nb_par);
			y:=[op(y),Y[1]];
			nb_par:=Y[2];  
		end if;
	end do;
	return([seq(y[i],i=1..nops(y))],time()-t0);
end proc:


#MPSupp_GRS(list,x,n,l,m)
#INPUT:  	list- [A_l(x),...,A_1(x),A_0(x)]- The list of the coefficients of the differential system. 
#		x- variable
#		n- Row dimension
#       	l- order of the differential system
#       	m- order of the regular solutions
#OUTPUT: 	the matrix polynomials [L_{1}(_mu),...,L_{m}(_mu)] associated to the differential system

MPSupp:=proc(list,x,n,l,m)
	local LL,i;
	for i from 1 to m do
		LL[i]:=add(ScalarMultiply(map(coeff,list[j],x,i),(_mu)^(l-j+1)),j=1..l+1);
	end do;
	return ([seq(LL[i], i = 1 .. m)]);
end proc:


#Second_Member_GRS(list,U,i,Maxk,lambda,n)
#INPUT:		list- The list of successive derivations of MPSupp- succdiffMPS
#		U- list of U_0, ..., U_{i-1}
#		i- The number of the coefficient U_i to compute
#		lambda- an eigenvalue
#		Maxk- Max{deg(U_j),j=0..i-1}+1 
#		n- Size of the differential system
#OUTPUT:	The list of the coefficients of -\sum_{j=0}^{i-1}L_{i-j}(Theta+lambda+j)U_{j} which is a polynomial in ln(x), i.e., [q_d, ..., q_0]

Second_Member:=proc(list,U,i,Maxk,lambda)
	local sm,o;
	for o from 1 to i do
		sm[o]:=[seq(add(simplify(Multiply(subs(_mu=lambda+o-1,list[i+1-o][k-j+1]),U[o][j])),j=1..k),k=1..Maxk)];
	end do;
	return([seq(simplify(-add(sm[o][j],o=1..i)),j=1..Maxk)]);
end proc:


#Inv_Solve_GRS(A,list1,b,j,U)
#INPUT:		A- MatrixInverse of L evaluated at _mu=lambda+i
#		list1- succdiff evaluated at _mu=lambda+i
#		b- the j-th element of Second_Member
#		U- A list of the computed coefficients of U_i
#OUTPUT:	Return the coefficient U_{i,p+1-j}

Inv_Solve_GRS:=proc(A,list1,b,j,U)
	local sm2,SM;
	sm2:=simplify(add(Multiply(list1[j-k],U[k]),k=1..j-1));
	SM:=b-sm2;
	return(simplify(Multiply(A,SM)));
end proc:


#Diff_non_integer_GRS(x,U0,max_kappa,ord,n,lambda,list1,list2,Inv)
#INPUT:		x- variable
#		lambda- the unique element of Ens[count] (an eigenvalue of L which does not differ by integers from the other eigenvalues)
#		kappa- the biggest partial multiplicity of lambda 
#		U0- The first polynomial vector in log(x) 
#		ord- The order of the power series involved in the regular solutions
#		n- Size of the differential system
#		list1- Succdiff
#		list2- SuccdiffMPS
#		Inv- Matrix inverse of L(_mu)
#OUTPUT: 	The general regular solution associated to the group Ens[count] or to lambda, where the power series is truncated up to the order ord.

Non_Integer_Difference:=proc(x,U0,ord,max_kappa,lambda,list1,list2,Inv)
	local U,i,j,k,sm,M,sd,seqU,UU,Y;
	seqU:=[U0];
	Y:=add(U0[i]*ln(x)^(max_kappa-i)/(max_kappa-i)!,i=1..max_kappa);
	for i from 1 to ord do
		sm:=Second_Member(list2,seqU,i,max_kappa,lambda);
		M:=convert(subs(_mu=lambda+i,Inv),Matrix);
		U[1]:=simplify(Multiply(M,sm[1]));
		UU:=U[1]*ln(x)^(max_kappa-1)/(max_kappa-1)!;
		sd:=[seq(subs(_mu=lambda+i,list1[k]),k=2..max_kappa)];
		for j from 2 to max_kappa do
			U[j]:=Inv_Solve_GRS(M,sd,sm[j],j,[seq(U[k],k=1..j-1)]);
			UU:=UU+U[j]*ln(x)^(max_kappa-j)/(max_kappa-j)!;
		end do;
		Y:=Y+UU*x^i;
		seqU:=[op(seqU),[seq(simplify(U[i]),i=1..max_kappa)]];
	end do;
	return(x^lambda*Y);
end proc:		



#LinearSolve_with_LU(A,b,n,set_parameters,nb_par)
#INPUT:		A- list of Matrices
#		b- List of Vectors
#		n- Size of each matrix
#		set_parameters- Parameters' set of the system
#OUTPUT:	Return the coefficient [U_{mp},..., U_{m0}] of U_m, when lambda+m is an eigenvalue of L(_mu)


LinearSolve_with_LU:=proc(A,b,n,set_parameters,nb_par)   
	local P,L,U,invP, Y,X,i,r,s,j,k,sol,tt,minus_indets, indets_b, ind, NP, proper_indets_Y;
	indets_b:=indets(b);
	NP:=nb_par;
	P,L,U:=LUDecomposition(A[1]);
	invP:=Transpose(P);
	r:=Rank(A[1]);
	Y[1]:=simplify(LinearSolve(L,invP.b[1],method='subs'));  # L is invertible
	X[1]:=simplify(LinearSolve(U,Y[1],method='subs',free=v));
	ind:=nops(indets(X[1]) minus `union`(indets_b,set_parameters));      # number of the constants v
	X[1]:=subs(seq(v[i]=c[NP+i],i=1..ind),X[1]);
	NP:=NP+ind;
	for i from 2 to nops(A) do
		Y[i]:=simplify(LinearSolve(L,invP.(b[i]-add(A[i-j].X[j+1],j=0..i-2)),method='subs'));
		proper_indets_Y:=indets(Y[i]) minus `union`(indets_b,set_parameters);
		sol:=solve({seq(Y[i][j],j=r+1..n)},proper_indets_Y);
		Y[i]:=simplify(subs(sol,Y[i]));
		for k from 1 to i-1 do
			X[k]:=simplify(subs(sol,X[k]));
		end do;
		#for j from r+1 to n do
		#	minus_indets:=indets(Y[i][j]) minus `union`(indets_b,set_parameters);
		#	if nops(minus_indets)<>0 then
		#		tt:=minus_indets[1];
		#		sol:=solve(Y[i][j],tt);
		#		Y[i]:=simplify(subs(tt=sol,Y[i]));
		#		for k from 1 to i-1 do
		#			X[k]:=simplify(subs(tt=sol,X[k]));
		#		end do;
		#	end if;
		#end do;
		X[i]:=simplify(LinearSolve(U,Y[i],method='subs',free=v));
		ind:=indets(X[i]) minus `union`(indets_b,set_parameters,{seq(c[i],i=1..NP)});
		X[i]:=subs(seq(ind[k]=c[NP+k],k=1..nops(ind) ),X[i]);
		NP:=NP+nops(ind);
	end do;
	return([seq(X[i],i=1..nops(A))],NP);
end proc:


#Diff_integer_GRS(x,U0,max_kappa,ord,n,lambda,list1,list2,list_diff,list_PM,Inv,set_parameters)
#INPUT: 	x- variable
#		lambda- the element of Ens[count] having the smallest real part 
#		max_kappa- the biggest partial multiplicity associated to lambda
#		U0- The first polynomial vector in log(x) 
#		ord- The order of the power series involved in the regular solutions
#		n- Size of the differential system
#		list1- Succdiff
#		list2- SuccdiffMPS
#		list_diff- A list containing the integer difference 
#		list_PM- A list containing the maximum values of the partial multiplicities of the eigenvalues differing from lambda by integers   
#		Inv- Matrix inverse of L
#		set_parameters- Parameters' set of the system
#OUTPUT: 	The general regular solution associated to the group Ens[count] where the power series is truncated up to the order max(list_diff)+ord.

Integer_Difference:=proc(x,U0,ord,max_kappa,n,p,lambda,list1,list2,list_diff,list_Mul,Inv,set_parameters,nb_par) 
	local sum_Mul, count, seqU, seqU_extend, i, sm, M, U, sd, j, Y, UU, sol, NP; 
	seqU:=[U0];
	NP:=nb_par;
	Y:=add(U0[i]*ln(x)^(max_kappa-i)/(max_kappa-i)!,i=1..max_kappa);	
	sum_Mul:=0;
	count:=0;
	seqU_extend:=[seqU[1]];
	for i from 1 to ord do
		if has(list_diff,i)=false then
			sm:=Second_Member(list2,seqU_extend,i,sum_Mul+max_kappa,lambda,n);
			M:=convert(subs(_mu=lambda+i,Inv),Matrix);
			U[1]:=simplify(Multiply(M,sm[1]));
			UU:=U[1]*ln(x)^(sum_Mul+max_kappa-1)/(sum_Mul+max_kappa-1)!;
			sd:=[seq(subs(_mu=lambda+i,list1[k]),k=2..sum_Mul+max_kappa)];
			for j from 2 to sum_Mul+max_kappa do
				U[j]:=Inv_Solve_GRS(M,sd,sm[j],j,[seq(U[k],k=1..j-1)]);
				UU:=UU+U[j]*ln(x)^(sum_Mul+max_kappa-j)/(sum_Mul+max_kappa-j)!;
			end do;
			Y:=Y+UU*x^i;
			seqU:=[op(seqU),[seq(U[k],k=1..sum_Mul+max_kappa)]];
			seqU_extend:=[op(seqU_extend),[seq(U[k],k=1..sum_Mul+max_kappa)]];
		else
			count:=count+1;
			sum_Mul:=sum_Mul+list_Mul[count];
			seqU_extend:=[seq([seq(Vector(n),kk=1..list_Mul[count]),op(seqU_extend[k])],k=1..i)];
			sm:=Second_Member(list2,seqU_extend,i,sum_Mul+max_kappa,lambda,n);
			sol:=LinearSolve_with_LU([seq(convert(subs(_mu=lambda+i,list1[k]),Matrix),k=1..sum_Mul+max_kappa)],sm,n,set_parameters,NP);
			U:=sol[1];
			NP:=sol[2];
			Y:=Y+add(U[k]*ln(x)^(sum_Mul+max_kappa-k)/(sum_Mul+max_kappa-k)!,k=1..nops(U))*x^i;
			seqU:=[op(seqU),U];
			seqU_extend:=[op(seqU_extend),U];		
		end if;
	end do;
	return(x^lambda*Y,NP);
end proc:

Verification:=proc(Liste_of_Matrices,Vector,x)
	local l, theta, V, list_thetaV, i;
	l:=nops(Liste_of_Matrices)-1;
	theta:=y->x*diff(y,x);
	list_thetaV:=[Vector];
	for i from 1 to l do
		list_thetaV:=[map(theta,list_thetaV[1]),op(list_thetaV)];
	end do;
	return(map(simplify,add(Liste_of_Matrices[i].list_thetaV[i],i=1..l+1)));
end proc: