#!/usr/bin/env python
# coding: utf-8

# # L'algorithme d'élimination de Gauss

# Commençons par écrire l'algorithme de Gauss standard. Il faut prendre garde au fait que les divisions qui apparaissent en cours de route nécessitent de travailler avec une matrice à coefficients dans un corps. Si ce n'est pas explicitement le cas au départ, on doit commencer par remplacer l'anneau R des coefficients de M, supposé intègre, par son corps des fractions K.

# In[1]:


def Gauss(M):
    N=copy(M)     #On renomme M afin de ne pas écraser cette matrice au cours de l'algorithme.
    L=N.nrows()
    C=N.ncols()
    R=M.base_ring()
    K=FractionField(R)
    N=N.change_ring(K) #On se place sur le corps des fractions de l'anneau de coordonnées de M.
    c=0
    while c<C:  #Travail sur la colonne c
        l=c
        while l<L and N[l,c]==0: #Recherche du premier coefficient non nul de cette colonne à partir de la ligne c.
            l=l+1
        if l==L:
            c=c+1
        else:
            if l!=c: 
                N.swap_rows(l,c) #Echange de lignes pour ramener le nouveau pivot en position (c,c).
            u=N[c,c]^(-1)    #Inversion du pivot
            for i in range(c+1,L):         
                N[i]=N[i]-N[i,c]*N[c]*u         #Elimination des coefficients sous le pivot.           
            c=c+1
    return N


# Un exemple numérique, qui met en évidence (sans surprise...) l'apparition de dénominateurs au cours de l'échelonnement.

# In[2]:


M=matrix(ZZ,[[4,1,-2,8],[5,-13,7,6],[-1,9,-11,9],[2,-8,12,-3]])
show(M), show(Gauss(M))


# Variante de l'algorithme d'élimination de Gauss sans division.

# In[3]:


def Gauss_entier(M):
    N=copy(M)
    L=N.nrows()
    C=N.ncols()
    c=0
    while c<C:
        l=c
        while l<L and N[l,c]==0:
            l=l+1
        if l==L:
            c=c+1
        else:
            if l!=c:
                N.swap_rows(l,c)
                
            for i in range(c+1,L):         
                N[i]=N[c,c]*N[i]-N[i,c]*N[c]     #Modification de la formule d'élimination pour supprimer les divisions.               
            c=c+1
    return N


# Un exemple numérique, qui suggère que les coefficients tendent à croître très fortement au cours de l'échelonnement.

# In[4]:


show(M), show(Gauss_entier(M))


# Un second exemple numérique, avec une matrice un peu plus grosse.

# In[5]:


MS=MatrixSpace(ZZ,8)
MM=MS.random_element()
show(MM)
show(Gauss_entier(MM))


# On peut également illustrer cette explosition de la taille des coefficients en travaillant dans ${\rm M}_n(K[X])$, où $K$ désigne un corps fini, et en observant la croissance du degré des coefficients dans l'algorithme.

# In[10]:


Z.<x>=PolynomialRing(ZZ)
M=random_matrix(Z,5); show(M)
show(Gauss_entier(M))


# # Décomposition LUP

# Réécrivons maintenant l'algorithme d'élimination de Gauss sans effectuer d'échange de ligne ; les pivots ne sont donc plus nécessairement diagonaux, et l'application de {1,...,n} dans lui-même associant à i l'indice de la ligne contenant le pivot de la colonne i est une permutation intrinsèquement associée à la matrice considérée.

# In[11]:


def Gauss_sans(M):
    N=copy(M)
    L=N.nrows()
    C=N.ncols()
    R=M.base_ring()
    K=FractionField(R)
    N=N.change_ring(K)
    c=0
    P=[]
    while c<C:
        l=0
        while l<L and (N[l,c]==0 or l in P): #Recherche du pivot dans la colonne c : il s'agit du premier 
            l=l+1                            #coefficient non nul (depuis le haut) ne figurant pas sur la ligne
        if l==L:                             #du pivot d'une colonne précédente.
            c=c+1
        else:
            P=P+[l]
            u=N[l,c]^(-1)    
            for i in range(l+1,L):         
                N[i]=N[i]-N[i,c]*N[l]*u      #Comme dans l'algorithme standard, on élimine les coefficients sous              
            c=c+1                            #le pivot.
    return N, P


# Exemple numérique :

# In[12]:


M=random_matrix(GF(7),6)
show(Gauss_sans(M))


# Calcul de la décomposition M=LUP

# In[13]:


def LUP(M):
    N=copy(M)
    C=N.ncols()
    p=[]
    R=N.base_ring()
    K=FractionField(R)
    N=N.change_ring(K)
    I=identity_matrix(K,C)
    L=copy(I)
    P=copy(I)
    c=0
    while c<C:
        l=0
        while l<C and (N[l,c]==0 or l in p):
            l=l+1
        if l==C:
            return "M n'est pas inversible"
        else:
            p=p+[l]
            u=N[l,c]^(-1)    
            for i in range(l+1,C):         
                L[:,l]=L[:,l]+N[i,c]*L[:,i]*u   #La matrice L est l'inverse du produit des matrices élémentaires
                N[i]=N[i]-N[i,c]*N[l]*u         #décrivant les opérations d'élimination successives. On l'obtient 
            c=c+1                               #en effectuant au fur et à mesure les opérations inverses sur les
                                                #colonnes de la matrice identité.
    for i in range(C):
        P[:,i]=I[:,p[i]]                        #La matrice de permutation P décrit la position des pivots
    
    U=N*P.transpose()
    return p,L,U,P, L*U*P==M


# Exemple numérique :

# In[17]:


M=random_matrix(GF(5),6)
show(LUP(M))


# Calcul de l'inverse d'une matrice L triangulaire inférieure unipotente, par échelonnement et report itéré des opérations effectuées sur la matrice identité. 

# In[18]:


def Inverse_TI(L):
    n=L.ncols()
    R=L.base_ring()
    K=FractionField(R)
    L=L.change_ring(K)
    I=identity_matrix(K,n)
    c=0
    while c<n:
        for i in range(c+1,n):
            I[i]=I[i]-L[i,c]*I[c] 
        c=c+1
    return I
            


# Exemple numérique :

# In[19]:


L=matrix(QQ,[[1,0,0,0],[2,1,0,0],[0,4,1,0],[2,2,0,1]])
show(L), show(Inverse_TI(L))


# Calcul de l'inverse d'une matrice triangulaire supérieure (non nécessairement unipotente).

# In[21]:


def Inverse_TS(U):
    n=L.ncols()
    R=U.base_ring()
    K=FractionField(R)
    U=U.change_ring(K)
    I=identity_matrix(K,n)
    c=n-1
    while c>-1:
        for i in range(0,c):
            u=U[c,c]^(-1)
            I[i]=I[i]-U[i,c]*I[c]*u
        I[c]=U[c,c]^(-1)*I[c]    
        c=c-1
    return I


# Exemple numérique :

# In[22]:


L=matrix(QQ,[[1,0,0,0],[2,1,0,0],[0,4,2,0],[2,2,0,1]])
U=L.transpose()
show(U),show(Inverse_TS(U))


# En utilisant conjointement la décomposition M=LUP et le calcul des inverses de L et U, on obtient immédiatement celui de l'inverse de M.

# In[28]:


def Inverse(M):
    n=M.nrows()
    R=M.base_ring()
    p,L,U,P,v=LUP(M)
    L1=Inverse_TI(L)
    U1=Inverse_TS(U)
    Q=P.transpose()
    N=Q*U1*L1
    return N, N*M==identity_matrix(R,n)


# Exemple numérique :

# In[29]:


M=random_matrix(QQ,4)
show(Inverse(M))


# Calcul du déterminant de M, à l'aide de la décomposition M=LUP.

# On commence par calculer ${\rm det}(P) = \varepsilon(p)$, où $p$ désigne la permutation dont $P$ est la matrice dans la base canonique de ${\bf K}^n$. Pour ce faire, on multiplie $p$ par des transpositions jusqu'à obtenir la permutation identité (c'est ainsi que l'on prouve que $p$ est un produit de transpositions).

# In[30]:


def signature(P):
    e=1
    n=P.nrows()
    c=0
    while c<n:
        i=c
        while P[i,c]==0:
            i=i+1
        if i==c:
            c=c+1
        else:
            P.swap_rows(0,i)
            e=-e
            c=c+1
            
    return  e   


# On peut alors calculer ${\rm det}(M)={\rm det}(U)\cdot{\rm det}(P)$.

# In[33]:


def determinant(M):
    p,L,U,P,v=LUP(M)
    n=M.nrows()
    u=1
    for i in range(n):
        u=u*U[i,i]
    e=signature(P)
    return e*u
    


# Un exemple numérique :

# In[34]:


M=random_matrix(ZZ,4)
show(M), show(determinant(M))


# # Exercice 1

# In[131]:


v=matrix([1,2,3,4]).transpose() ; show(v)


# In[132]:


M=matrix([[1,2,3],[4,5,6]]); show(M)


# In[133]:


M.parent() 


# Voici une matrice "aléatoire" de taille 5x5 à coefficients entiers

# In[134]:


M=random_matrix(ZZ,5); show(M) 


# Les exemples suivants montrent des matrices 10x10 dont les coefficients sont donnés par différentes variables aléatoires iid : loi uniforme sur (0,1), loi gaussienne centrée réduite, loi uniforme sur {-10,-9,...,9,10} et loi uniforme sur un corps fini à 9 éléments.

# In[135]:


L=[random() for i in range(100)] 
M=matrix(10,L); show(M)    


# In[136]:


L=[gauss(0,1) for i in range(100)]
M=matrix(10,L) ; show(M)


# In[137]:


L=[randint(-10,10) for i in range(100)] 
M=matrix(10,L) ; show(M)


# In[158]:


K.<a>=GF(9)
K.modulus()


# Remarque : si q est une puissance d'un nombre premier p, la commande GF(q) produit un corps fini à q éléments. On peut obtenir le polynôme minimal du générateur a par la commande K.modulus() ou charpoly(a) (polynôme caractéristique de la multiplication par a dans GF(q) vu comme un espace vectoriel sur GF(p)).

# In[159]:


L=[randint(0,2)+randint(0,2)*a for i in range(100)]
M=matrix(10,L);show(M)


# Et voici enfin une matrice 3x3 dont les coefficients sont des indéterminées.

# In[162]:


R.<T11,T12,T13,T21,T22,T23,T31,T32,T33>=ZZ[]
M=matrix(3,[T11,T12,T13,T21,T22,T23,T31,T32,T33]) ; show(M),show(det(M))


# La commande MatrixSpace(A,n,m) définit l'espace des matrices de taille mxn à coefficients dans l'anneau A.

# In[165]:


MS=MatrixSpace(ZZ,4)
M=MS.random_element(); show(M)


# In[166]:


I=MS.one(); show(I)


# In[167]:


charpoly(M)


# In[168]:


minpoly(M)


# In[169]:


MS.base_ring()


# In[170]:


N=M.change_ring(GF(5)); show(N)


# In[171]:


N.base_ring()


# Opérations algébriques usuelles sur les matrices : multiplication par un vecteur, multiplication matricielle, inversion et transposition.

# In[174]:


show(N*v)


# In[176]:


show(N^4), show(N^(-1)), show(N.transpose())


# # Exercice 2

# In[178]:


A=matrix(QQ,[[1,2,-1,2], [1,1,1,1],[-1,2,1,2],[-1,2,1,0]]);show(A)
I=identity_matrix(QQ,4); show(I)
B=block_matrix([[I,A],[I,I]]) ; show(B)
C=block_matrix([[I,I],[I,I]]) ; show(C)


# In[ ]:


B.right_kernel()


# In[180]:


B.column_space()


# In[181]:


C.right_kernel()


# In[182]:


C.column_space()


# In[183]:


B.echelon_form()


# In[184]:


L=matrix(QQ,[[3,4,5]])


# In[187]:


L.right_kernel()


# In[188]:


LL=matrix([[3,4,5]])
LL.right_kernel()


# In[189]:


D,U,V = LL.smith_form()
show(D,U,V)


# # Exercice 3

# Soit $A \in {\rm GL}_n({\bf C})$. En appliquant le procédé d'orthonormalisation de Gram-Schmidt aux colonnes de $A$, nous obtenons une matrice triangulaire supérieure $R$ à coefficients diagonaux dans ${\bf R}_{>0}$ telle que les colonnes de $Q=AR$ forment une base unitaire de ${\bf C}^n$ : la matrice $Q$ est donc unitaire, et nous obtenons ainsi la décomposition $A=QR^{-1}$. On en déduit l'identité $$|\det(A)| = \det(R^{-1}) = \det(R)^{-1}.$$
# Par ailleurs, le coefficient diagonal $R_{i,i}$ est l'inverse de la norme du vecteur $V_i$ déduit de la $i$-ème colonne $A(i)$ de $A$ par orthogonalisation relativement aux colonnes précédentes ; on a $A(i) = V_i + W_i$, avec $W_i \perp V_i$, donc $$||A(i)||^2 = ||V_i||^2 + ||W_i||^2 \geqslant ||V_i||^2,$$ c'est-à-dire $$R_{i,i}^{-1} \leqslant ||A(i)||$$
# et $$|\det(A)| = |{\rm det}(R^{-1})| = \prod_{i=1}^n R_{i,i}^{-1} \leqslant \prod_{i=1}^n ||A(i)||.$$

# Calcul de la borne de Hadamard d'une matrice carrée à coefficients complexes :

# In[35]:


def Borne_Hadamard(M):
    M.change_ring(CC) #On se place sur le corps des complexes
    n=M.nrows()
    N=M.transpose()
    b=1
    for i in range(n):
        b=b*norm(N[i])   #Si v est un vecteur dans C^n, norm(v) est le carré de sa norme hermitienne standard.
    return ceil(b.n())   #Pour faciliter l'utilisation de la borne de Hadamard, on se contente de calculer 
                         #sa partie entière supérieure. 
    
    


# Un exemple numérique:

# In[36]:


M=random_matrix(ZZ,3)
show(M);
show(Borne_Hadamard(M));


# On peut illustrer l'intérêt de la borne de Hadamard par le \emph{calcul modulaire} des déterminants. Commençons par observer qu'il est naturel de s'attendre à ce que le déterminant d'une matrice carrée $M$ de taille $n \times n$ \`a coefficients entiers soit en général de l'ordre de grandeur de $B^n$, où B désigne la maximum des valeurs absolues des coefficients de M : c'est en effet la somme de $n!$ termes qui sont chacun des produits de $n$ coefficients de $M$. Illustrons-le numériquement, avec une matrice $10 \times 10$ dont les coefficients sont compris entre $-10000$ et $1000$ ; son déterminant doit avoir environ $\lfloor {\rm log}_{10}\left((10^4)^{10}\right)\rfloor = 40$ chiffres.

# In[48]:


L=[randint(-10^4,10^4) for i in range(100)] 
M=matrix(10,L) ; show(M)
d=det(M) ; show(d)
floor(numerical_approx(log(abs(d),10)))


# Il est possible de calculer un tel déterminant en conservant des matrices dont tous les coefficients restent de taille raisonnable, disons comparable à celle des coefficients de $M$ : c'est le $calcul$ $modulaire$. Il suffit pour cela de choisir des nombres premiers $p_1,\ldots, p_r$ de cette taille, de calculer le $r$ déterminants $d_i={\rm det}(M \hskip1mm ({\rm mod}~p_i)$, identifiés à des entiers entre $0$ et $p_i-1$, puis finalement d'invoquer le th\'eor\`eme chinois des restes pour reconstituer ${\rm det}(M)$. Si $P = p_1 \cdots p_r$ est tel que $P>2H$, où $H$ désigne la borne de Hadamard de $M$, alors $d={\rm det}(M)$ est l'unique nombre entier compris entre $-\frac{P}{2}$ et $\frac{P}{2}$ tel que $d \equiv d_i \hskip1mm ({\rm mod}~p_i)$ pour tout $i=1,\ldots, r$.

# In[57]:


P=Primes()
L=[P[i] for i in range(1600,1611)]; show(L)
numerical_approx(log(prod(L),10))
show(prod(L));
show(floor(numerical_approx(log(prod(L),10))));
prod(L)>2*Borne_Hadamard(M)


# In[59]:


D=[det(M.change_ring((GF(L[i])))) for i in range(0,11)]; show(D)
E=[d.lift() for d in D] ;
show(crt(E,L)-prod(L));


# # Exercice 4

# In[60]:


R.<X11,X12,X13,X21,X22,X23,X31,X32,X33>=ZZ[]
M=matrix(R,3,[X11,X12,X13,X21,X22,X23,X31,X32,X33]); show(M)


# In[61]:


M[1]=X11*M[1]-X21*M[0]
M[2]=X11*M[2]-X31*M[0]
show(M)


# In[62]:


M[2]=M[1,1]*M[2]-M[2,1]*M[1]
show(M)


# In[63]:


show(factor(M[2,2]))


# In[64]:


def Gauss_Bareiss_Z(M):
    N=copy(M)
    L=N.nrows()
    C=N.ncols()
    c=0
    d=1
    while c<C:
        l=c
        while l<L and N[l,c]==0:
            l=l+1
        if l==L:
            c=c+1
        else:
            if l!=c:
                N.swap_rows(l,c)
                
            for i in range(c+1,L):  
                for j in range(c+1,C):
                    N[i,j]=(N[c,c]*N[i,j]-N[i,c]*N[c,j])//d
                N[i,c]=0   
            d=N[c,c]   
            c=c+1
    return N,d


# In[65]:


M=matrix(ZZ,[[4,1,-2,8,1],[5,-13,7,6,1],[-1,9,-11,9,1],[2,-8,12,-3,1]])
show(Gauss_entier(M)), show(Gauss_Bareiss_Z(M))


# In[66]:


def Gauss_Bareiss(M):
    N=copy(M)
    L=N.nrows()
    C=N.ncols()
    R=M.base_ring()
    K=R.fraction_field()
    c=0
    d=1
    while c<C:
        l=c
        while l<L and N[l,c]==0:
            l=l+1
        if l==L:
            c=c+1
        else:
            if l!=c:
                N.swap_rows(l,c)
                
            for i in range(c+1,L):  
                for j in range(c+1,C):
                    N[i,j]=(N[c,c]*N[i,j]-N[i,c]*N[c,j])//d
                N[i,c]=0   
            d=N[c,c]   
            c=c+1
    return N,d


# In[67]:


R.<X>=GF(7)[]
A=random_matrix(GF(7),6)
M=X*identity_matrix(GF(7),6)-A ; show(M)


# In[68]:


show(MM);
show(Gauss_Bareiss(MM))