» Home Page
Topologia Algebrica (2011-2S)
--» Forma Normale di Smith, MCD e Omologia (2011-2S)

Classe Matrice (matrice.py)

La prima cosa da fare è scrivere una classe per rappresentare in modo sintetico le operazioni sulle matrici. Un esempio è il seguente.

class Matrice:
 def __init__(self,data,uno=1):
  """Inizializzazione di una matrice, in cui le righe
vengono completate da zeri se necessario.
Il parametro opzionale uno=1.0 oppure uno=1
e' per avere interi, oppure float, oppure...
Esempio: Matrice([[1,2],[1]])
"""
  if type(data) != type([]):
   data=[[data]]
  self.ismatrix=True
  self.data=data
  self.uno=uno
  self.num_righe=len(data)
  self.num_colonne=max([len(r) for r in data])
  for d in self.data:
   if len(d) < self.num_colonne:
    d += ([0] * (self.num_colonne-len(d)))
 def __repr__(self):
  res="[\n"+("\n".join([str(riga) for riga in self.data]))+"\n]"
  return str(res)
 
 def __str__(self):
  return  self.__repr__()
 
 def __getitem__(self,key):
   "per avere A[i,j]"
   return self.data[key[0]][key[1]]

 def __setitem__(self,key,value):
   "quando A[i,j] = 2"
   self.data[key[0]][key[1]]=value
 
 def colonna(self,j):
  return [ self[i,j] for i in range(self.num_righe) ]

 def riga(self,i):
  return [ self[i,j] for j in range(self.num_colonne) ]
 
 def sottomatrice(self,ii,jj):
  """sottomatrice (ii,jj)-esimo"""
  return Matrice ( [ [ self.data[i][j] 
                    for j in range(self.num_colonne) if j != jj
                  ] for i in range(self.num_righe) if i != ii ])

 def __or__(self,other):
  """ A|B : matrice a blocchi aumentata """
  if self.num_righe != other.num_righe:
   raise Exception ("numero righe non coincide")
  return Matrice ( 
  [self.data[i] + other.data[i] for i in range(self.num_righe)])

 def __div__(self,other):
  """ A/B: matrice a blocchi uno sopra l'altro"""
  if self.num_colonne != other.num_colonne:
   raise Exception("numero di colonne non coincide")
  return Matrice (self.data + other.data)
 
 def t(self):
  """trasposta"""
  return Matrice ( [ [ self.data[i][j] 
                    for i in range(self.num_righe)
                  ] for j in range(self.num_colonne) ])

 def __mod__(self,other):
  """coefficienti mod other"""
  return Matrice( [[ (self.data[i][j] % other ) 
        for j in range (self.num_colonne)
      ] for i in range(self.num_righe) ])
                 


 def __add__(self,other):
  return Matrice ( [ [ 
        (self.data[i][j] + other.data[i][j]) 
        for j in range (self.num_colonne)
      ] for i in range(self.num_righe) ])

 def __eq__(self,other):
  for i in range(self.num_righe): 
   for j in range (self.num_colonne):
       if (self.data[i][j] != other.data[i][j]):
        return False
  return True

 def __sub__(self,other):
  """A-B"""
  return Matrice ([ [ (self.data[i][j] - other.data[i][j]) 
                       for j in range (self.num_colonne)
                     ] for i in range(self.num_righe) ])

 def __mul__(self,other):
  """ self * other : prodotto di matrici, righe per colonne """
  return Matrice ( 
                  [ [ sum ([self.data[i][k] * other.data[k][j] 
                     for k in range(self.num_colonne) 
                          ])  for j in range(other.num_colonne)
                   ] for i in range(self.num_righe) ])

 def __rmul__(self,other):
  """numero * matrice """
  return Matrice( 
                 [ [ self.data[i][j] * other 
                    for j in range(self.num_colonne)
                  ] for i in range(self.num_righe) ])

 def __neg__(self):
  """- matrice """
  return Matrice( 
               [ [ - self.data[i][j] 
                   for j in range(self.num_colonne)
                 ] for i in range(self.num_righe) ])


 def inversa(self):
  """inversa"""
  if self.num_colonne==1:
   return self.uno / self[1,1]
  else:
   return self.uno/det(self) * \
    trasposta( 
     Matrice( [ [(-1)**(i+j) * det(self.sottomatrice(i,j))  
                 for j in range(self.num_colonne) 
               ] for i in range(self.num_righe) ] ))
 

 def __pow__(self,other):
  """ A ** n potenza. """
  result= MatriceIdentica(self.num_righe)
  if other >= 0:
   for i in range(other):
    result = self * result
   return result
  else:
   invself=self.inversa()
   for i in range(-other):
    result = invself  * result
   return result 
 
 def __xor__(self,other):
    """ A ^ n : potenza (sinonimo di A^n ) """
    print "attenzione alle precedenze su ^!"
    return self.__pow__(other)


def MatriceIdentica(n):
 x=MatriceNulla(n,n)
 for i in range(x.num_colonne):
  x.data[i][i] = 1
 return x

def MatriceNulla(n,m):
 return Matrice([[0 for j in range(m)] for i in range(n)])

def MatriceDiagonale(l):
 le=len(l)
 x=MatriceNulla(le,le)
 for i in range(le):
  x.data[i][i] = l[i]
 return x
 
def trasposta(m):
  return m.t()

def det(x):
  """determinante di x"""
  if x.num_colonne!=x.num_righe:
   raise Exception("Matrice non quadrata")
  if x.num_colonne==1: return x.data[0][0]
  return sum ([(-1)**i * x.data[0][i]*det(x.sottomatrice(0,i))  
              for i in range(x.num_colonne)  ]) 

def Estendi(m):
 return ( Matrice([ [1] ]) | MatriceNulla(1,m.num_colonne)
        ) / ( MatriceNulla(m.num_righe,1) | m )


def MatricePermutazione( nupla ):
 n=len(nupla)
 m=MatriceNulla(n,n)
 m.uno=1
 for i in range(n): m[i,nupla[i]] = 1
 return m  

def Pij(i,j,n):
 p=range(n)
 p[i]=j
 p[j]=i
 return MatricePermutazione(p)

Calcolo di MCD (mcd.py)

La funzione MCD calcola il massimo comun divisore di una riga $a=(a_1,\ldots, a_n)$ di interi, insieme a una matrice $R$ tale che \[ aR = (MCD,0,0,\ldots, 0) \] Questo ci consente (con la prima colonna di $R$) di risolvere la equazione di Bezout \[ k_1a_1 + \ldots + k_n a_n = MCD \]

Per esempio, per risolvere il problema dei marinai e delle noci di cocco occorre risolvere l'equazione diofantea \[ 1024 n - 15625 k = 11529, \] dove $n$ è il numero di noci di cocco iniziali.

from mcd import *

a=[1024,-15625]
l,R=MCD(a)
print l
print R
==>

[-1L, 0L]
[
[4776L, 15625L]
[313L, 1024L]
]

Quindi \[ 1024 \times 4776 - 15625 \times 313 = -1, \] da cui segue che una soluzione dell'equazione è \[ (n,k) = ( 4776 \times (-11529) , 313 \times (-11529) ), \] e tutte le soluzioni si ottengono sommando multipli della sconda colonna della matrice $R$ trovata sopra, per cui tutte le soluzioni sono \[ n \in \{ - 4776 \times 11529 + h 15625 : h \in \mathbb{Z} \}. \] La soluzione con $n$ positivo minimo è \[ n = -4776 \times 11529 \mod 15625 = 15621 \]

Il codice è questo:

from matrice import *

def EuclidMatrix(b):
 n=len(b) + 1
 M=MatriceIdentica(n)
 for i in range(1,n):
  M[0,i] = b[i-1] 
 return M

def srcMCD(a,R):
 n = len(a) 
 if len( [x for x in a[1:]  if x != 0] ) == 0:
   # se il primo e' l'unico nullo: termina.
   return (a,R)
 imin=0
 valmin=max( [abs(x) for x in a] )
 for i in reversed(range(n)):
  if (a[i] != 0 ) and (abs(a[i]) <= valmin) :
    imin = i
    valmin = abs(a[i])
 permM=Pij(0,imin,n)
 a[0],a[imin] = a[imin],a[0]
 if len( [x for x in a[1:]  if x != 0] ) == 0 :
  return ( a , R* permM )
 euclM = EuclidMatrix([ - a[i]/a[0] for i in  range(1,n) ] )
 c = Matrice([a]) * euclM
 return srcMCD(c.data[0], R*permM * euclM )

def MCD(a):
 a=[1L*x for x in a]
 n=len(a)
 return srcMCD(a,MatriceIdentica(n))


if __name__=='__main__':
 tmpa=[12*450*123, 15*3*123, 454*5*123, 35*5*12*123,50*12*9*123] 
 l,R=MCD(tmpa)
 print tmpa, "*", R , "=", Matrice([tmpa]) * (R)
 print "=?= ", l


Forma normale di Smith (smith.py)

Con un po' di sforzo possiamo utilizzare la funzione MCD per determinare la forma normale di Smith di una matrice di interi. Si tratta di un algoritmo dimostrativo, incredibilmente lento ma di cui è facile dimostrare che termina sempre.

from matrice import *
from mcd import *
import random
VERBOSE=False

def MinAij(M):
 """ posizione del coefficiente minimo nonzero """
 ijmin=[0,0]
 valmin=max( max([abs(x) for x in M.riga(j)]) for j in range(M.num_righe) ) 
 for i in reversed(range(M.num_righe)):
  for j in reversed(range(M.num_colonne)):
   if (M[i,j] != 0 ) and (abs(M[i,j]) <= valmin) :
    ijmin = [i,j]
    valmin = abs(M[i,j])
 return ijmin


def srcSmith(L,M,R):
 if VERBOSE: print "chiamato: ", M
 n=M.num_righe
 m=M.num_colonne
 # first find the smallest non-zero 
 i,j = MinAij(M)
 val=M[i,j]
 if i != 0 or j != 0:
  PL = Pij(0,i,n)
  PR = Pij(0,j,m)
  M = PL * M * PR
  L=PL * L
  R = R * PR
  if VERBOSE: print "permuted: ", M
 a=M.riga(0)
 if m>1 and len( [x for x in a[1:]  if x != 0] ) != 0 :
   ap,Ra = MCD(a)
   M = M * Ra
   R = R * Ra
   if VERBOSE: print "after a", M
 b=M.colonna(0)
 if n>1 and len( [x for x in b[1:]  if x != 0] ) != 0 :
   bp,Rb = MCD(b)
   #trasposta
   M = Rb.t() * M
   L = Rb.t() * L
   if VERBOSE: print "after b", M
 if n==1 or m == 1:
  return (L,M,R)
 MatMod = M % M[0,0]
 if MatMod  == MatriceNulla(n,m) and (len( [x for x in b[1:]  if x != 0] ) == 0 ) and ( len( [x for x in a[1:]  if x != 0] ) == 0 ) :
  xL,xM,xR = Smith(M.sottomatrice(0,0))
  return (
   Estendi(xL) * L, 
   ( Matrice([M.riga(0)]) / ( MatriceNulla(n-1,1) | xM ) ), 
   R * Estendi(xR) 
         )
 elif MatMod == MatriceNulla(n,m):
  if VERBOSE: print "matmod == 0 "
  return srcSmith(L,M,R)
 else:
  ti,tj = MinAij(MatMod)
  for i in range(n):  
   M[i,0] += M[i,tj]
  for j in range(m):  
   R[j,0] += R[j,tj]
  return srcSmith(L,M,R)
 return Fail

def Smith(M):
 n=M.num_righe
 m=M.num_colonne
 L,SM,R = srcSmith(MatriceIdentica(n),M,MatriceIdentica(m))
 # sistemiamo i segni
 mm=MatriceIdentica(m)
 for i in range(min(n,m)):
  if SM[i,i] < 0 :
   SM[i,i] = - SM[i,i]
   mm[i,i] = -  mm.uno
 return L, SM, R * mm

def test(NMAX):
 n=random.randint(20,40)
 m=random.randint(20,40)
 M=Matrice( [[random.randint(-NMAX,NMAX) *1L for i in xrange(m)] for j in xrange(n) ])
 print "nxm=", n, "x", m
 L,SM,R = Smith(M)
 return (L * M * R == SM ) 

if __name__=='__main__':
 M=Matrice([[6L,0L,9L],[0L,5L,6L]])
 M=Matrice([[6,0,9,14,15,16],[0,5,6,12,13,14]])
 M=Matrice([[3,2,3],[0,2,0],[2,2,2] ])
 L,SM,R= Smith(M)
 print L,SM,R
 for x in range(100):
  print "x= ", x
  xx=test(100)
  if not xx:
   raise Exception("failed!\n")
 M=Matrice([[1024 , -3125] ])
 L,SM,R= Smith(M)
 print SM


Vero Software

Esiste in letteratura software già pronto e ottimizzato per calcolare MCD e forma normale di Smith. Per esempio, con pari/gp, basta usare la funzione matsnf per la forma normale (con l'ordinamento decrescente dei divisori elementari non nulli, cioè dei termini sulla diagonale della matrice in forma normale).

? M=[3,2,3;0,2,0;2,2,2]
%1 = 
[3 2 3]

[0 2 0]

[2 2 2]

? SM=matsnf(M,flag=1)
%2 = [[2, 1, -3; 0, 0, 1; 1, 0, 0], [-1, -2, 1; 0, 3, -1; 1, 0, 0], [0, 0, 0; 0, 2, 0; 0, 0, 1]]

Omologia (hom.py)

Avendo a disposizione un qualsiasi algoritmo per calcolare la forma normale di Smith, è finalmente possibile calcolare l'omologia di un complesso di catene (sugli interi) in modo molto semplice:

from matrice import *
from smith import *


def Hom(cpl):
 """cpl e' la lista delle matrici dei bordi (cosa succede se uno dei C_k e' nullo?"""
 N=len(cpl)+1
 cpl=[Matrice([[ 1 for x in range(cpl[0].num_righe)]])]+cpl
 allR = [ None for x in range(N) ]
 allL = [ None for x in range(N) ]
 allSM =[ None for x in range(N) ]
 allSM[0] = Matrice([ [0 for x in range(cpl[0].num_colonne) ] ]) 
 for k in range(1,N):
  allL[k], allSM[k], allR[k] = Smith(cpl[k])
 p = [ [allSM[k][x,x] for x in \
  range(min(allSM[k].num_righe,allSM[k].num_colonne)) \
  if allSM[k][x,x] != 0] for k in range(N) ]
 p += [ [] ] # k = N+1
 c = [ allSM[k].num_colonne for k in range(N) ]
 return [ [ c[k]-len(p[k]) - len(p[k+1]) , [d for d in p[k+1] if d != 1] ] \
 for k in range(N) ]

def pretty(l):
 f=l[0]
 d=l[1] 
 r=[]
 if f>0:
  r += [ "Z^%i" % f ]
 if len(d)>0:
  r += [ "Z_%i" % x for x in d ]
 if r == [] :
  r = [ "0" ]
 return "+".join(r) 

if __name__=='__main__':
 cpl = [ Matrice( [[0]]), Matrice([[2]]) ]  #proiettivo
 cpl = [ Matrice( [[0,0]]), Matrice([[2],[0]]) ]  #klein
 cpl = [ Matrice( [[0,0]]) , Matrice ([[0],[0]]) ] #toro 
 hom=Hom(cpl)
 N=len(hom)
 for k in range(N):
  print  "H_%i=" %k , pretty(hom[k])
   

Ovviamente gli algoritmi Poincaré-Smith possono essere estesi ad anelli diversi da $\mathbb{Z}$, e ottimizzati per il calcolo dell'omologia (cfr. e.g. [http://www.sagemath.org/doc/reference/sage/homology/chain_complex.html]).