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

# # Exercice 1

# In[1]:


v=vector([1,2,3,4]) ; show(v)


# In[2]:


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


# In[3]:


M=matrix(10,lambda i,j : i+j+1) ; show(M)


# In[4]:


M.parent() 


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

# In[5]:


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[6]:


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


# In[7]:


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


# In[8]:


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


# In[8]:


K.<a>=GF(9)
L=[randint(0,2)+randint(0,2)*a for i in range(100)]
M=matrix(10,L);show(M)


# 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[9]:


M.parent()


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

# In[10]:


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(R,n,m) définit l'espace des matrices de taille mxn à coefficients dans l'anneau R.

# In[24]:


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


# In[25]:


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


# In[20]:


MS.base_ring()


# In[14]:


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


# In[15]:


N.base_ring()


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

# In[16]:


show(N*v)


# In[17]:


show(v*N)


# In[18]:


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


# # 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[27]:


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[29]:


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


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

# In[31]:


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[35]:


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


# Soit dit en passant, Sage contient déjà une commande permettant d'échelonner une matrice donnée. Appliquée à notre matrice à coefficients entiers, il semble que l'algorithme utilisé permette d'éviter les dénominateurs tout en conservant des coefficients de taille raisonnable (cf. exercice 4)...

# In[37]:


show(M.echelon_form())


# # Exercice 2

# 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[41]:


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[44]:


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


# Calcul de la décomposition M=LUP

# In[45]:


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)
    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
    P=I.matrix_from_columns(p)                  #colonnes de la matrice identité.
    U=N*P.transpose()
    return L,U,P, L*U*P==M


# Exemple numérique :

# In[49]:


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


# In[26]:





# 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[50]:


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[53]:


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))


# In[ ]:


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


# In[86]:


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[91]:


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[95]:


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


# Exemple numérique :

# In[99]:


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[157]:


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[159]:


def determinant(M):
    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[167]:


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


# # 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[178]:


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[177]:


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


# # Exercice 4

# In[1]:


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[2]:


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


# In[3]:


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


# In[4]:


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


# In[5]:


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[6]:


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


# In[ ]: