############################################################################

# This file contains a list of procedures related to my ISSAC'03 paper   
# "Factorization of Differential Systems in Characteristic $p$ and Application". 

# Needed Maple packages are:
# - linalg, 
# - DEtools, 
# - tensor (only for the "eigmodp" file).

#----------------------------------------------------------------------------
 
# List of the procedures:

# LiePol(A,n,p): compute the $p$-curvature (and all elements of the Lie sequence) 
# of the system Y'=A Y of dimension n with coefficients in F_p(z)

# ---

# % Char_Pol_Modp(M,lambda,z,p): compute the characteristic polynomial \chi(lambda) 
# of a matrix M with entries in F_p(z).

# ---

# % LBTM([A1,...Ar],r,n) (resp. UBTM): form the Lower (resp. Upper) Block triangular 
# Toeplitz Matrix of size nr with given square blocks A[i] of size n.

# ---

# % pdecomposition(f,p): decompose an element f of F_p(z) on the basis (1,z,...,z^{p-1}) 
# of F_p(z) over F_p(z^p): f = f_0 + f_1 z + ... + f_{p-1} z^{p-1}, f_i in F_p(z^p). 

# ---

# % pdecomposition_mat(A,n,p): decompose the matrix A of size n on the basis (1,z,...,z^{p-1}) 
# of F_p(z) over F_p(z^p): A = A_0 + A_1 z + ... + A_{p-1} z^{p-1}, A_i matrix of size n with 
# coefficients in F_p(z^p).

# ---
 
# % diff_mat(A,p): compute the derivative w.r.t. z of a matrix A with coeffs in F_p(z)

# ---

# % matrice_op(A,n,p): compute the matrix of the operator D = d/dz - A in the basis form by the   
# Eij where Eij= vect( 0$(i-1), z^j , 0$(n-i) )(i=1..n,j=0..p-1).

# ---

# % matrice_op_bp(A,n,p): compute the matrix N (see section 1.3 of the paper) of d/dz - A in 
# a particular basis.    

# ---

# % rat_sol1(A,n,p): compute the rational solutions of the system Y'=A Y where A is a matrix of size
# n with coefficients in F_p(z) by computing the Nullspace of matrice_op(A,n,p).

# ---

# % rat_sol(A,n,p): compute the rational solutions of the system Y'=A Y where A is a matrix of size
# n with coefficients in F_p(z) by computing the Nullspace of matrice_op_bp(A,n,p).

# ---

# % AbsFact(f,p,mult): factor the polynomial f with coefficients in F_p(z^p) 
# over \overline{F}_p(z^p).

# ---

# % ScalpCurv(A,n,p): diagonalize the system Y'=A Y (where A is a matrix of size n with coefficients
# in F_p(z)) having a scalar $p$-curavture. (see section 2.1.4 of the paper) 

# ---

# % IsoDec(A,Pcurv,p,opt): find the isotypical decomposition of the system Y' =A Y having 
# $p$-curvature Pcurv. We compute the characteristic spaces of Pcurv. (see section 2.3 of the paper)
# The decomposition is computed over F_p if opt=0 and over \overline{F}_p if opt=1.

# ---

# % ScalDec(A,n,p,opt): compute the scalar decomposition of Y'=A Y. We first compute 
# IsoDec(A,Pcurv,p,opt) and then we diagonalize blocks with scalar $p$-curvature (see section 2.4 
# of the paper). The option opt has the same goal as in the IsoDec procedure.

# ---

# % eigmodp(A,n,p) compute a basis (over F_p(z^p) of the eigenring of Y'=AY where A is a matrix 
# of size n with coefficients in F_p(z). 

# ------------------------------------------------------------------------------

# For all details, see my ISSAC'03 paper "Factorization of Differential
# Systems in Characteristic $p$ and Application". 

# For all questions, write to "thomas.cluzeau@unilim.fr".


###############################################################################

LiePol := proc(A,n,p)
 
# INPUT: the matrix A (size n) of the differential system mod p
# OUTPUT: the p-curvature (a matrix of size n)
 
# Compute the $p$-curvature of a differential system 
# by computing one by one the elements of the associated 
# Lie sequence without denominators. This is an improvment of Katz' method in 
# the sense that we make matrices multiplications faster by doing denominator free computations.
 
local i, j, B, S, De, d1;

# we use global variables for elements of the Lie sequence and the denominator of A
# because we need them to compute a basis of rational solutions with Katz' formula 
# in the procedure ScalpCurv

global _Lie_Pol, _d:
 
De:=map(X->denom(X),A); 
_d:=lcm(seq(seq(De[i,j],j=1..n),i=1..n)): # the denominator of the matrix

d1:=diff(_d,z):

B:=evalm(_d*A): # the matrix with polynomial coefficients

_Lie_Pol[0]:=Matrix(n,n,shape=identity):

_Lie_Pol[1]:=-B:
 
for i from 2 to p do

# we adapt the classic formula A_{i+1} = A_i ' - A A_i to have denominator free computations

       _Lie_Pol[i]:=map(X->modp(Normal(X),p),evalm(evalm(_d*diff_mat(_Lie_Pol[i-1],p))- 
B &* _Lie_Pol[i-1]-evalm((i-1)*d1*_Lie_Pol[i-1])));

 od:

 RETURN(map(X->modp(Normal(X),p),evalm(_Lie_Pol[p]/(_d^p)))):

 end:


########################################################################

Char_Pol_modp := proc(M,lambda,z,p)

# Input: matrix M with coefficients in F_p(z), the variable of the output lambda, the
# variable z and the prime p 
# Output: the characteristic polynomial of the matrix M (a polynomial in F_p[z,lambda] 
# where the constant z^p is replaced by c) 

# This procedure applies in particular when M is the $p$-curavture and in this case
# the output is in F_p[c][lambda] with c:=z^p.

# Compute the characteristic polynomial of a matrix which entries are rational functions 
# in two variables with coefficients in the finite field F_p using a denominator free 
# interpolation method. The "existing" Det(M-lambda I) mod p" procedure is not efficient 
# because of the two variables. Consequently we use an interpolation method to reduce the 
# problem to computing determinants of matrices with polynomials entries in only one 
# variable. 

# I thank M. van Hoeij for the idea to improve this characteristic polynomial computation 
# with an interpolation approach 

# WARNING : in the input, p must be greater than n

local n, Id, Mat, i, S, X, Mat_pol:

n:=rowdim(M):

Id:=Matrix(n,n,shape=identity):

# we form the characteristic matrix Mat

Mat:=evalm(lambda*Id-M):

# we eliminate the denominators columns by columns

Mat_pol := map( X -> simplify( X ) mod p , concat ( seq ( evalm ( simplify ( col ( Mat , i ) * 
lcm ( seq ( denom ( col ( Mat , i ) [ j ] ) , j = 1..n ) ) ) ) , i = 1..n ) ) ) ; 

X:=[seq(i,i=0..n)]:

# we specialize the variable \lambda= 0,..,n and we compute determinants of 
# matrices with polynomial in only one variable as coefficients. 
# Note that we need n+1 distinct interpolation points which implies the restriction p > n.

S:= [seq( Det ( evalm ( map ( v -> subs ( lambda = X[i] , v ) , Mat_pol ) ) ) mod p , 
i= 1..n+1)]:

# we reconstruct by interpolation with the existing Interp() mod p procedure

RETURN( collect(Rem(Primpart(Interp ( X , S , lambda ) mod p,lambda) mod p,z^p-c,z) mod p,
lambda)):

end:

##############################################################################

LBTM:=proc(L,r,n,p)

# Input: L =[A1,.. Ar ] a list of  nxn matrices
# Output:  nrxnr  lower triangular Toeplitz matrix 

local N,j,T;
if r > 0 then 
N:=matrix(n,n, 0);
T:=augment(stackmatrix(op(L)),seq(stackmatrix(N$j,op(L)[1..r-j]),j=1..r-1));
fi;
end: 

###############################################################################

UBTM:=proc(L,r,n)

# Input: L =[A1,.. Ar ] a list of  nxn matrices
# Output:  nrxnr  lower triangular Toeplitz matrix 

local N,j,T;
if r > 0 then 
N:=matrix(n,n, 0);
T:=stackmatrix(augment(op(L)),seq(augment(N$j,op(L)[1..r-j]),j=1..r-1));
fi;
end: 

################################################################################

pdecomposition := proc(f,p)

# Input: an element f of F_p(z)
# Output: a sequence of p elements of F_p(z^p)

# Computes the decomposition of f on the
# basis 1,..,z^(p-1) of F_p(z) over F_p(z^p)

local PP,i,f1,Q,P1,S1,R,Q1,j:

# we separate the numerator and the denominator
PP:= numer(f):
Q := denom(f):

# we replace the constant z^p by c

P1:=Rem(PP,z^p-c,z) mod p: 
Q1:=Rem(Q,z^p-c,z) mod p:

# we compute the inverse S1 of the denominator Q1 using Bezout method

Gcdex(Q1,z^p-c,z,'S1','t') mod p:	 

# so that our element becomes the following R

R:=Rem(P1*S1,z^p-c,z) mod p:

# we extract each coefficient

for j from 0 to (p-1) do
	f1[j+1]:= coeff(R,z,j) mod p:
od:

RETURN(seq(f1[i],i=1..p)):

end:

###################################################################################

pdecomposition_mat := proc(A,n,p)

# In : an n*n matrix A with coeffs in F_p(z)
# Out: a sequence of p n*n matrices with coeffs in F_p(z^p)

# Compute the decomposition of a matrix A on the
# basis 1,..,z^(p-1) of F_p(z) over F_p(z^p)
# we only aplly the above procedure to each elements of A

local CC,i,Mat:

CC:=convert(map(pdecomposition,A,p,alpha),array):
for i from 1 to p do 
	Mat[i]:=concat(seq(col(CC,i+j*p),j=0..(n-1))): 
od:

RETURN(seq(evalm(Mat[i]),i=1..p)):

end:

######################################################################

diff_mat := proc(A,p)

# In : a square matrix A with coeffs in F_p(z) 
# Out: a square matrix  with coeffs in F_p(z) (same size as A)

# Compute the derivative d/dz of the matrix A  

RETURN(map(modp,map(diff,evalm(A),z),p));     

end:

#######################################################################

matrice_op := proc(A,n,p)

# In : an n*n matrix A with coeffs in F_p(z)
# Out: an (n*p)*(n*p) matrix with coeffs in F_p(z^p) 

# Compute the matrix of the operator D = d/dz - A in the basis of the   
# Eij where Eij= vect( 0$(i-1), z^j , 0$(n-i) )(i=1..n,j=0..p-1).
  
local i,j,k,L,l,d,E,M,w,v,v1,v2,R1,q,r:

M:=matrix(p*n,p*n):
v2:=vector(n):

for i from 1 to p*n do
    d:=iquo(i,p,'r'):
    if r <> 0 then
       q:=d+1:
    else
       q:=d:	
    fi:         	   
    v:=vector(n,0):
    r:=(i-1) mod p:
    v[q]:=z^r:       
    v2:=evalm(map(proc(f) diff(f,z) end proc,v)-multiply(A,v)):
    L[i]:=convert(map(pdecomposition,v2,p,alpha),array);
od:

M:=concat(seq(L[i],i=1..p*n)):

RETURN(M):

end:

########################################################################

mat_op_bp := proc(A,n,p)

# In : an n*n matrix A with coeffs in F_p(z)
# Out: an (n*p)*(n*p) matrix with coeffs in F_p(z^p) 

# Compute the matrix of the operator D = d/dz - A
# in a particular basis where its shape seems interesting.

# It consists in decomposing A and the unknown rational solution 
# in the basis 1,...,z^(p-1)

# Remark: the matrix N of our "paper"

local SS, TildeN, Z, Id, L, Mat_Der, N, i:

SS:=[pdecomposition_mat(A,n,p)];

TildeN:=evalm(LBTM([seq(SS[i],i=1..p)],p,n)+UBTM([diag(0$n),
seq(c*SS[p-i],i=0..p-2)],p,n)):

Z:=Matrix(n,n,shape=zero):
Id:=Matrix(n,n,shape=identity):

for i from 1 to p-1 do 
	L[i]:=concat(seq(Z,j=1..i),-i*Id,seq(Z,j=i+2..p)):
od:
L[p]:=concat(seq(Z,i=1..p)):
Mat_Der:=map(modp,stackmatrix(seq(L[i],i=1..p)),p):

N:=map(modp,evalm(TildeN+Mat_Der),p):

RETURN(N):

end:

####################################################################

rat_sol1 := proc(A,n,p)

# In : an n*n matrix A with coeffs in F_p(z)
# Out: the sequence containing the rational solutions
# of the system Y'=AY and { } if there are not.  

# Compute the rational solutions of the system Y'=AY
# by solving the linear system of size n*p given by the matrix 
# operator D = d/dz - A in the basis of the E_ij (see matrice_op)

local B, N1, N, N2, i, j, l, k, De, d:

B := evalm(matrice_op(A,n,p)):
De:=map(X->denom(X),B):
d:=lcm(seq(seq(De[i,j],j=1..n*p),i=1..n*p)):
N1:=Nullspace(evalm(d*B)) mod p:
N2:=matrix(nops(N1),n):
k:=0:
for i from 1 to nops(N1) do
   for j from 1 to n do
   N2[i,j] := add(N1[i][l+k+1]*z^l,l=0..(p-1)):
   k:=k+p:
   od:
   N[i]:=[seq(N2[i,j],j=1..n)]:
   k:=0:
od:

RETURN(seq(N[i],i=1..nops(N1)));   

end:

################################################################

rat_sol2 := proc(A,n,p)

# In : an n*n matrix A with coeffs in F_p(z)
# Out: the sequence containing the rational solutions
# of the system Y'=AY and { } if there are not.  

# Compute the rational solutions of the system Y'=AY
# by solving the linear system of size n*p given by the matrix 
# operator D = d/dz - A in an interesting particular basis 
# (see matrice_op_bp)

local i, j, NN, N1, s, y, De, d:

NN:=evalm(mat_op_bp(A,n,p)):
De:=map(X->denom(X),NN):
d:=lcm(seq(seq(De[i,j],j=1..n*p),i=1..n*p)):
N1:=Nullspace(evalm(d*NN)) mod p:

for i from 1 to nops(N1) do
   for j from 0 to p-1 do
   	y[i,j] := vector(n,[seq(N1[i][k+j*n],k=1..n)]);
   od:
od:

for i from 1 to nops(N1) do
	s[i]:=add(y[i,j]*z^j,j=0..p-1):
od:   

RETURN(seq(evalm(s[i]),i=1..nops(N1)));   

end:

##############################################################

AbsFact := proc(f,p,mult)

# In : a polynomial f in F_p(z)[T] irreducible over F_p(z)
# Out: a sequence of factors in \bar{F_p}(z)[T] given like 
# that :[factor, multiplicity]

# The method used consists in computing first the polynomial that will define
# (by a RootOf) the extension of F_p needed and then trying to factor f absolutely 
# in this extension using the "Factor(f,RootOf) mod p" command. 

local test, n, m, g, S, minim, R, i:

test:=0:
n:=degree(f,T):
m:=degree(f,z):
g:=gcd(n,m):

# If gcd(degree(f,z),degree(f,T))=1 then the polynomial is absolutely irreducible 

if g=1 then

   RETURN([f,mult]);

else

# else we define the polynomial that will define (by a RootOf) the extension
	while (test=0) do
	      minim:=randpoly(x,degree=g,coeffs=rand(0..p-1)):
	    
	      if ((Irreduc(minim) mod p)=true and degree(minim,x)=g and lcoeff(minim) = 1) then
		 test:=1:
	      else
		 test:=0:	 
	      fi:
	od: 

# we try to factor f in the suitable extension

	print(minim,f);
	S:=Factors(f,RootOf(minim)) mod p:	
fi:

for i from 1 to nops(op(2,S)) do
	R[i]:=[op(1,op(i,op(2,S))),mult*op(2,op(i,op(2,S)))]:
od:
	

RETURN(seq(R[i],i=1..nops(op(2,S))));

end:

#######################################################################################

ScalpCurv := proc(A,n,p)

# Input: a square matrix A such that A_{p} = - \lambda I, the size n of A and the prime p
# Output: the diagonal matrix B = \mu I such that [B] is equivalent to [A]

local i, mu, Id, B, Inter:
global _PP, _IPP, _M:

mu := simplify(-trace(A)/n) mod p:  # satisfies \tau(mu) = lambda

Id:=Matrix(n,n,shape=identity):

B:=evalm(A+mu*Id): # B_p = A_p + \tau(mu) = 0

LiePol(B,n,p):

# we compute a fundamental matrix  _PP of rational solutions of B with Katz' method

_PP := map(X->Normal(X) mod p,evalm(add((-z)^i/factorial(i)*_Lie_Pol[i]/_d^i,i=0..p-1))):
_IPP := Inverse(_PP) mod p:

# we return _PP[A] = [-mu I]

Inter:=map(X->modp(Normal(X),p),evalm(_IPP &* A &* _PP)):
_M:=map(X->modp(Normal(X),p),evalm(Inter - _IPP &* diff_mat(_PP,p))):

RETURN(_M):

end:

#################################################################################

IsoDec := proc(A,Pcurv,p,opt)

# In : a square matrix A of size n, its pcurvature Pcurv, the prime p and an option 
# opt which is 0 or 1
# Out: a matrix of size n i.e "the isotypical decomposition" of Y'=AY

# We keep the following elements as global variables as this procedure is destined to
# be a "subprocedure" of our procedure ScalDec for the "maximal" factorisation of our 
# system:

# _Nbre1 := number of factors in the isotypical decomposition
# _Nbre2 := the sequence of the dimensions of the factors
# _PP := the matrix of passage (_IPP := its inverse)
# _MMM := The system in its form after the isotypical factorisation
#  _G1 := the sequence of the distincts absolutly irreducible factors of the 
# characteristic polynomial of Pcurv with their multiplicity

# if opt = 0 then we factor over F_p(z) and if opt=1 we factor over \overline{F}_p(z)

local F, i, G, N, inter, j, l, GG, GGG:
global _Nbre1, _Nbre2, _PP, _MMM, _G1, _IPP:

# We compute and factor over F_p(z) the characteristic polynomial of Pcurv

F:=op(2,Factors(Char_Pol_modp(Pcurv,T,z,p)) mod p):

GGG:=select(X->has(X,T),F):

# If opt=1 and we have one factor of degree greater than 1, we try an absolute 
# factorisation of it

for i from 1 to nops(GGG) do
		
	if (degree(GGG[i][1],T) > 1) and opt = 1 then
		G[i]:=AbsFact(GGG[i][1],p,GGG[i][2]):	
	else
		G[i]:=GGG[i]:
	fi:
od:

_G1:=[seq(G[i],i=1..nops(GGG))]:

_Nbre1:=nops(_G1):

# _Nbre1 = 1 implies that the system does not factor if we only look at the characteristic 
# polynomial of the pcurvature

if _Nbre1=1 then 
_MMM:=A:
_PP:=Matrix(rowdim(A),shape=identity):
_IPP:=_PP:
_Nbre2[1]:=rowdim(A):

else

# _Nbre1 > 1 means that we have at least two factors; we thus compute the associated characteristic spaces
# _Nbre2[i] = size of the ist factor = dimension of the associated characteristic space

for i from 1 to _Nbre1 do
	N[i]:=Nullspace(evalm(subs(T=Pcurv,subs(c=z^p,_G1[i][1]^_G1[i][2])))) mod p:
	_Nbre2[i]:=nops(N[i]):
od:

# We deduce the associated matrix of passage 

_PP:=concat(seq(seq(convert(N[i][j],matrix),j=1..nops(N[i])),i=1.._Nbre1)):
_IPP:=Inverse(_PP) mod p:

# we compute then the factorisation of the system

inter:=map(X->modp(Normal(X),p),evalm(_IPP &* A &* _PP)):
_MMM:=map(X->modp(Normal(X),p),evalm(inter - _IPP &* diff_mat(_PP,p))):

fi:

RETURN(_MMM):

end:
				

#####################################################################################

ScalDec := proc(A,n,p,opt)

# In : an n*n matrix A with entries in F_p(z)
# Out: an n*n matrix with entries in \bar{F_p}(z)

# We factor the system Y'=AY by using only the informations given by the pcurvature

local k,i,Pcurv,Id,Z,Id1,c,so,so1,c1,Pcurvbis,Pcurv1,B,A1,Z1,G2,multip,m:

Id:=Matrix(n,n,shape=identity):
Z:=Matrix(n,n,shape=zero):

# we compute the pcurvature of the system

Pcurv := LiePol(A,n,p):

# c:=Pcurv[1,1]:

# We test if it is scalar (c*Id); if it is we solve the equation
# y^p+y^(p-1) = -c and we have our factorisation. 

if IsMatrixShape(convert(map(X->modp(Normal(X),p),Pcurv),Matrix),scalar) then
	
 		RETURN((Normal(-trace(A)/n) mod p)*Id):
	
else	

# if it is not, we decompose the system using the characteristic spaces of its pcurvature
# for this we apply the "dec_int" procedure
 	
	evalm(IsoDec(A,Pcurv,p,opt)):

# If _Nbre1=1, then the characteristic polynomial of the pcurvature
# is some F^m with m>=1. If m=1 then we return A which is irreducible. 
# it is our irreducibility criterion

	if ((_Nbre1=1) and (op(2,op(1,_G1))=1)) then 
		print(irreducible):
		RETURN(A):
	else

# Else we have a non trivial isotypical decomposition of the system
# and we deduce the pcurvature associated with each factor

	Pcurvbis:=map(X->modp(Normal(X),p),evalm(_IPP &* Pcurv &* _PP)):
	k:=1:

	for i from 1 to _Nbre1 do
		A1[i]:=convert(_MMM,Matrix)[k..k+_Nbre2[i]-1,k..k+_Nbre2[i]-1]:
		Pcurv1[i]:=map(X->modp(Normal(X),p),convert(Pcurvbis,Matrix)[k..k+_Nbre2[i]-1,k..k+_Nbre2[i]-1]):
		k:=k+_Nbre2[i]:	
	od:

# we try now to decompose each factor
	for i from 1 to _Nbre1 do
	
# _Nbre2[i]=1 means that the factor is of dimension 1 thus irreducible	
		if _Nbre2[i]=1 then
			B[i]:=evalm(A1[i]):
		else

# else we test whether the pcurvature is scalar or not; if it is we "redecompose" the factor
# into a diagonal one 

			
			Id1[i]:=Matrix(rowdim(Pcurv1[i]),rowdim(Pcurv1[i]),shape=identity):
			if IsMatrixShape(convert(map(X->modp(Normal(X),p),Pcurv1[i]),Matrix),scalar) then
			
					B[i]:=evalm((Normal(-trace(A1[i])/_Nbre2[i]) mod p)*Id1[i]):	
			
			else
                            
				B[i]:=evalm(A1[i]):
			fi:
			
		fi:
	
		
	od:

fi:	
		
		
		
	RETURN(evalm(diag(seq(map(X->modp(Normal(X),p),evalm(B[i])),i=1.._Nbre1)))):
			
fi:

end:	

####################################################################################################################

eigmodp := proc(A,n,p)

# In : an n*n matrix A with coeffs in F_p(z)
# Out: a sequence of n*n matrices that constitutes a basis of the
# eigenring of the system Y'=AY

# Remark : this method needs to solve a linear system of size n^2*p !! 

local i,j,k,l,TA,TRANSPA,ID,TENS,r,MAT,N,MAT1,RATSOL,Outp,Ligne:

r:=1:

# we "tensorise" to obtain a system of the form Z'=MAT1 Z
# where MAT1 = A°I-I°transp(A) (° stands for the tensor product of matrices)

TA := create([1,1],convert(A,array)):
TRANSPA :=create([1,1],convert(transpose(A),array)):
ID:=create([1,1],convert(Matrix(1..n,1..n,shape=identity),array)):
TENS:=lin_com(prod(TA,ID),-1,prod(ID,TRANSPA)):
if op(2,op(1,op(2,op(1,TENS)))) = [1,1,1,1]  then
  MAT:=op(2,op(2,op(2,op(1,TENS)))):
else
  MAT:=op(2,op(1,op(2,op(1,TENS)))):
fi:  
N:=array(1..n^2):
for i from 1 to n do
    for j from 1 to n do
        N[r]:=[seq(seq(MAT[i,k,j,l],l=1..n),k=1..n)]:
        r:=r+1:
   od:
od:             

MAT1:=map(X->modp(Normal(X),p),evalm(stackmatrix(seq(N[i],i=1..n^2))));

# we compute the matrix of the linear operator d/dz - MAT1
# in order to find the rational solutions of the system Y'= MAT1 Y

RATSOL := [rat_sol2(MAT1,n^2,p)]:
Outp:=array(1..nops(RATSOL)):
for i from 1 to nops(RATSOL) do
        
	k:=1:
	for j from 1 to n do

		Ligne[j] := [seq(RATSOL[i][j],j=k..n+k-1)]:
	        k:=k + n: 

	od:
	Outp[i]:=stackmatrix(seq(Ligne[j],j=1..n)):
od:

RETURN(seq(map(X->modp(Normal(X),p),evalm(Outp[i])),i=1..nops(RATSOL)));

end:

########################################################################################################################