"Fossies" - the Fresh Open Source Software Archive

Member "NZMATH-1.2.0/nzmath/matrix.py" (19 Nov 2012, 64187 Bytes) of package /linux/misc/old/NZMATH-1.2.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "matrix.py" see the Fossies "Dox" file reference documentation.

    1 from __future__ import division
    2 
    3 import nzmath.ring as ring
    4 import nzmath.vector as vector
    5 
    6 
    7 class Matrix(object):
    8     """
    9     Matrix is a class for matrices.
   10     """
   11 
   12     def __init__(self, row, column, compo=0, coeff_ring=0):
   13         """
   14         Matrix(row, column [,components, coeff_ring])
   15         """
   16         self._initialize(row, column, compo, coeff_ring)
   17         self._selectMatrix()
   18 
   19     def _initialize(self, row, column, compo=0, coeff_ring=0):
   20         """
   21         initialize matrix.
   22         """
   23         if (isinstance(row, (int, long))
   24             and isinstance(column, (int, long))
   25             and row > 0
   26             and column > 0 ): # row and column check
   27             self.row = row
   28             self.column = column
   29             self.compo = []
   30             if isinstance(compo, (ring.Ring, int)):
   31                 coeff_ring = compo
   32                 compo = 0
   33             if compo == 0:
   34                 zero = 0
   35                 if coeff_ring != 0:
   36                     coeff_ring = _getRing(coeff_ring)
   37                     zero = coeff_ring.zero
   38                 for i in range(self.row):
   39                     self.compo.append([zero] * self.column)
   40             else:
   41                 if isinstance(compo[0], list):
   42                     if len(compo) != self.row:
   43                         raise ValueError("len(compo) " + 
   44                                          "is not match the matrix size")
   45                     for i in range(self.row):
   46                         if len(compo[i]) != self.column:
   47                             raise ValueError("len(compo[%d]) " % i + 
   48                                              "is not match the matrix size")
   49                         self.compo.append(compo[i][:])
   50                 elif isinstance(compo[0], tuple):
   51                     if len(compo) != self.column:
   52                         raise ValueError("len(compo) " +
   53                                          "is not match the matrix size")
   54                     self.compo = [[] for i in range(self.row)]
   55                     for i in range(self.column):
   56                         if len(compo[i]) != self.row:
   57                             raise ValueError("len(compo[%d]) " % i + 
   58                                              "is not match the matrix size")
   59                         j = 0
   60                         for ele in compo[i]:
   61                             self.compo[j].append(ele)
   62                             j += 1
   63                 elif isinstance(compo[0], vector.Vector):
   64                     if len(compo) != self.column:
   65                         raise ValueError("len(compo) " +
   66                                          "is not match the matrix size")
   67                     self.compo = [[] for i in range(self.row)]
   68                     for i in range(self.column):
   69                         if len(compo[i]) != self.row:
   70                             raise ValueError("len(compo[%d]) " % i +
   71                                              "is not match the matrix size")
   72                         j = 0
   73                         for ele in compo[i].compo:
   74                             self.compo[j].append(ele)
   75                             j += 1
   76                 else:
   77                     if (len(compo) != self.row * self.column):
   78                         raise ValueError("len(compo) " +
   79                                          "is not match the matrix size")
   80                     for i in range(self.row):
   81                         self.compo.append(
   82                         compo[self.column * i : self.column * (i + 1)])
   83             if coeff_ring == 0:
   84                 self.coeff_ring = ring.getRing(self.compo[0][0])
   85             else:
   86                 self.coeff_ring = coeff_ring = _getRing(coeff_ring)
   87                 if not(isinstance(ring.getRing(self.compo[0][0]), coeff_ring.__class__)):
   88                     compos = []
   89                     for i in range(self.row):
   90                         compos.append(map(coeff_ring.createElement, self.compo[i]))
   91                     self.compo = compos
   92         else:
   93             raise ValueError("invalid value for matrix size")
   94 
   95     def _selectMatrix(self):
   96         """
   97         Select Matrix class.
   98         """
   99         if self.coeff_ring.isfield():
  100             if self.row == self.column:
  101                 self.__class__ = FieldSquareMatrix
  102             else:
  103                 self.__class__ = FieldMatrix
  104         else:
  105             if self.row == self.column:
  106                 self.__class__ = RingSquareMatrix
  107             else:
  108                 self.__class__ = RingMatrix
  109 
  110     def __getitem__(self, index):
  111         """
  112         M[i,j] : Return (i,j)-component of M.
  113         M[i] <==> M.getColumn(i)
  114         """
  115         if isinstance(index, tuple):
  116             return self.compo[index[0] - 1][index[1] - 1]
  117         elif isinstance(index, (int, long)):
  118             return self.getColumn(index)
  119         else:
  120             raise IndexError("Matrix invalid index: %s" % index)
  121 
  122     def __setitem__(self, key, value):
  123         """
  124         M[i,j] = a   :   Substitute a for (i,j)-component of M.
  125         M[i] = V <==> M.setColumn(i, V)
  126         """
  127         if isinstance(key, tuple):
  128             self.compo[key[0] - 1][key[1] - 1] = value
  129         elif isinstance(key, (int, long)):
  130             self.setColumn(key, value)
  131         else:
  132             raise TypeError(self.__setitem__.__doc__)
  133 
  134     def __eq__(self, other):
  135         """
  136         Check self == other.
  137         self == 0 means whether self == zeromatrix or not.
  138         """
  139         if isinstance(other, Matrix):
  140             if (self.row != other.row) or (self.column != other.column):
  141                 return False
  142             for i in range(self.row):
  143                 for j in range(self.column):
  144                     if self.compo[i][j] != other.compo[i][j]:
  145                         return False
  146             return True
  147         elif isinstance(other, int) and other == 0: # zero matrix ?
  148             return not bool(self)
  149         else:
  150             return False
  151 
  152     def __hash__(self):
  153         val = 0
  154         for i in range(self.row):
  155             for j in range(self.column):
  156                 val += hash(self.compo[i][j])
  157         return val
  158 
  159     def __ne__(self, other):
  160         """
  161         Check self != other.
  162         self != 0 means whether self != zeromatrix or not.
  163         """
  164         return not (self == other)
  165 
  166     def __nonzero__(self):
  167         """
  168         Check self != zeromatrix
  169         """
  170         for i in range(self.row):
  171             for j in range(self.column):
  172                 if self.compo[i][j]:
  173                     return True
  174         return False
  175 
  176     def __contains__(self, item):
  177         """
  178         Check whether item == one of self's element.
  179         """
  180         for i in range(self.row):
  181             if item in self.compo[i]:
  182                 return True
  183         return False
  184 
  185     def __repr__(self):
  186         return_str = ""
  187         for i in range(self.row):
  188             return_str += str(self.compo[i])
  189             if i + 1 != self.row:
  190                 return_str += "+"
  191         return return_str
  192 
  193     def __str__(self):
  194         return_str = ""
  195         width = [1] * self.column      # width of each column
  196         for j in range(self.column):
  197             for i in range(self.row):
  198                 if len(str(self.compo[i][j])) > width[j]:
  199                     width[j] = len(str(self.compo[i][j]))
  200         for i in range(self.row):
  201             for j in range(self.column):
  202                 return_str += "%*s " % (width[j], str(self.compo[i][j]))
  203             return_str += "\n"
  204         return return_str[:-1]
  205 
  206     def __call__(self, arg):
  207         """
  208         Return matrix applied __call__ to each elements.
  209         """
  210         sol = []
  211         for i in range(self.row):
  212             for j in range(self.column):
  213                 ele = self.compo[i][j]
  214                 if callable(ele):
  215                     sol.append(ele(arg))
  216                 else:
  217                     sol.append(ele)
  218         return createMatrix(self.row, self.column, sol)
  219 
  220 
  221     def map(self, function):
  222         """
  223         Return matrix applied function to all self elements.
  224         """
  225         compos = []
  226         for i in range(self.row):
  227             compos = compos + (map(function, self.compo[i]))
  228         return createMatrix(self.row, self.column, compos)
  229 
  230     def reduce(self, function, initializer=None):
  231         """
  232         Return value applied binomial function to all self elements.
  233         """
  234         if initializer:
  235             val = reduce(function, self.compo[0], initializer)
  236         else:
  237             val = reduce(function, self.compo[0])
  238         for i in range(1, self.row):
  239             val = reduce(function, self.compo[i], val)
  240         return val
  241 
  242     def copy(self):
  243         """
  244         Create a copy of the instance.
  245         """
  246         compos = []
  247         for i in range(self.row):
  248             for j in range(self.column):
  249                 compos.append(self.compo[i][j])
  250         Mat = self.__class__(self.row, self.column, compos, self.coeff_ring)
  251         return Mat
  252     
  253     def set(self, compo):
  254         """
  255         set(compo) : Substitute list for components
  256         """
  257         if (len(compo) != self.row * self.column):
  258             raise ValueError("len(compo) is not match the matrix size")
  259         for i in range(self.row):
  260             self.compo[i] = compo[self.column*i : self.column*(i + 1)]
  261 
  262     def setRow(self, m, arg):
  263         """
  264         setRow(m, new_row) : new_row should be a list/Vector
  265         """
  266         if isinstance(arg, list):
  267             if (len(arg) != self.column):
  268                 raise vector.VectorSizeError(
  269                     "len(compo) is not match the row size")
  270             self.compo[m - 1] = arg[:]
  271         elif isinstance(arg, vector.Vector):
  272             if (len(arg) != self.column):
  273                 raise vector.VectorSizeError(
  274                     "len(compo) is not match the row size")
  275             self.compo[m - 1] = arg.compo[:]
  276         else:
  277             raise TypeError(self.setRow.__doc__)
  278 
  279     def setColumn(self, n, arg):
  280         """
  281         setColumn(n, new_column) : new_column should be a list/Vector
  282         """
  283         if isinstance(arg, list):
  284             if (len(arg) != self.row):
  285                 raise ValueError("len(compo) is not match the column size")
  286             for i in range(self.row):
  287                 self.compo[i][n - 1] = arg[i]
  288         elif isinstance(arg, vector.Vector):
  289             if (len(arg) != self.row):
  290                 raise ValueError("len(compo) is not match the column size")
  291             for i in range(self.row):
  292                 self.compo[i][n - 1] = arg.compo[i]
  293         else:
  294             raise TypeError(self.setColumn.__doc__)
  295 
  296     def getRow(self, i):
  297         """
  298         getRow(i) : Return i-th row in form of Matrix
  299         """
  300         return vector.Vector(self.compo[i - 1][:])
  301 
  302     def getColumn(self, j):
  303         """
  304         getColumn(j) : Return j-th column in form of Matrix
  305         """
  306         column = []
  307         for k in range(self.row):
  308             column.append(self.compo[k][j - 1])
  309         return vector.Vector(column)
  310 
  311     def swapRow(self, m1, m2):
  312         """
  313         swapRow(m1, m2) : Swap self's m1-th row and m2-th row.
  314         """
  315         tmp = self.compo[m1 - 1][:]
  316         self.compo[m1 - 1] = self.compo[m2 - 1][:]
  317         self.compo[m2 - 1] = tmp[:]
  318 
  319     def swapColumn(self, n1, n2):
  320         """
  321         swapColumn(n1, n2) : Swap self's n1-th column and n2-th column.
  322         """
  323         for k in range(self.row):
  324             tmp = self.compo[k][n1 - 1]
  325             self.compo[k][n1 - 1] = self.compo[k][n2 - 1]
  326             self.compo[k][n2 - 1] = tmp
  327 
  328     def insertRow(self, i, arg):
  329         """
  330         insertRow(i, new_row) : added new_row
  331         new_row can be a list or a Matrix
  332         """
  333         if isinstance(arg, list):
  334             if self.column != len(arg):
  335                 raise vector.VectorSizeError()
  336             self.compo.insert(i - 1, arg)
  337             self.row += 1
  338         elif isinstance(arg, vector.Vector):
  339             if self.column != len(arg):
  340                 raise vector.VectorSizeError()
  341             self.compo.insert(i - 1, arg.compo)
  342             self.row += 1
  343         elif isinstance(arg, Matrix):
  344             if self.column != arg.column:
  345                 raise MatrixSizeError()
  346             self.compo += arg.compo
  347             self.row += arg.row
  348         else:
  349             raise TypeError()
  350         self._selectMatrix()
  351 
  352     def insertColumn(self, j, arg):
  353         """
  354         insertColumn(j, new_column) : added new_column
  355         new_column can be a list or a Matrix
  356         """
  357         if isinstance(arg, list):
  358             if self.row != len(arg):
  359                 raise vector.VectorSizeError()
  360             for k in range(self.row):
  361                 self.compo[k].insert(j-1, arg[k])
  362             self.column += 1
  363         elif isinstance(arg, vector.Vector):
  364             if self.row != len(arg):
  365                 raise vector.VectorSizeError()
  366             for k in range(self.row):
  367                 self.compo[k].insert(j-1, arg.compo[k])
  368             self.column += 1
  369         elif isinstance(arg, Matrix):
  370             if self.row != arg.row:
  371                 raise MatrixSizeError()
  372             for k in range(arg.row):
  373                 self.compo[k] = \
  374                 self.compo[k][:j - 1] + arg.compo[k] + self.compo[k][j - 1:]
  375             self.column += arg.column
  376         else:
  377             raise TypeError()
  378         self._selectMatrix()
  379 
  380     def extendRow(self, arg):
  381         """
  382         extendRow(new_row) : join new_row in vertical way.
  383         """
  384         self.insertRow(self.row + 1, arg)
  385 
  386     def extendColumn(self, arg):
  387         """
  388         extendColumn(new_column) : join new_column in horizontal way.
  389         """
  390         self.insertColumn(self.column + 1, arg)
  391 
  392     def deleteRow(self, i):
  393         """
  394         deleteRow(i) : deleted i-th row
  395         """
  396         self.row -= 1
  397         del self.compo[i - 1]
  398         self._selectMatrix()
  399 
  400     def deleteColumn(self, j):
  401         """
  402         deleteColumn(j) : deleted j-th column
  403         """
  404         self.column -= 1
  405         for k in range(self.row):
  406             del self.compo[k][j - 1]
  407         self._selectMatrix()
  408 
  409 
  410     def transpose(self):
  411         """
  412         Return transposed matrix of self.
  413         """
  414         trans = []
  415         for j in range(1, self.column + 1):
  416             for i in range(1, self.row + 1):
  417                 trans.append(self[i, j])
  418         return self.__class__(self.column, self.row, trans, self.coeff_ring)
  419 
  420     def getBlock(self, i, j, row, column=None):
  421         """
  422         Return a block whose size is row*column, (1,1)-element is self[i,j].
  423         """
  424         if column == None:
  425             column = row
  426         if i + row - 1 > self.row or j + column - 1 > self.column:
  427             raise MatrixSizeError()
  428         mat = []
  429         for k in range(i - 1, i + row - 1):
  430             mat.extend(self.compo[k][j - 1:j + column - 1])
  431         return createMatrix(row, column, mat, self.coeff_ring)
  432 
  433     def subMatrix(self, I, J=None):
  434         """
  435         Return submatrix whose element is self[i, j] for i in I and j in J.
  436         If I, J is not index(list or tuple) but integer,
  437          return submatrix which deleted I-th row and J-th column from self.
  438         """
  439         if J == None:
  440             J = I
  441         if isinstance(I, (int, long)) and isinstance(J, (int, long)):
  442             mat = self.copy()
  443             mat.deleteRow(I)
  444             mat.deleteColumn(J)
  445             return mat
  446         else:
  447             mat = []
  448             for i in I:
  449                 for j in J:
  450                     mat.append(self[i, j])
  451             return createMatrix(len(I), len(J), mat, self.coeff_ring)
  452 
  453     def toMatrix(self, flag=True):
  454         """
  455         return the matrix (self)
  456         (this method is for compatibility to vector)
  457         """
  458         return self
  459 
  460 
  461 class SquareMatrix(Matrix):
  462     """
  463     SquareMatrix is a class for square matrices.
  464     """
  465 
  466     def __init__(self, row, column=0, compo=0, coeff_ring=0):
  467         """
  468         SquareMatrix(row [, column ,components, coeff_ring])
  469         SquareMatrix must be row == column .
  470         """
  471         self._initialize(row, column, compo, coeff_ring)
  472         self._selectMatrix()
  473 
  474     def _initialize(self, row, column=0, compo=0, coeff_ring=0):
  475         """
  476         initialize matrix.
  477         """
  478         if isinstance(compo, (ring.Ring, int)):
  479             coeff_ring = compo
  480             compo = 0
  481         if isinstance(column, list):
  482             compo = column
  483             column = row
  484         elif isinstance(column, ring.Ring):
  485             coeff_ring = column
  486             column = row
  487         elif column == 0:
  488             column = row
  489         if row != column:
  490             raise ValueError("not square matrix")
  491         Matrix._initialize(self, row, column, compo, coeff_ring)
  492 
  493     def isUpperTriangularMatrix(self):
  494         """
  495         Check whether self is upper triangular matrix or not.
  496         """
  497         for j in range(1, self.column + 1):
  498             for i in range(j + 1, self.row + 1):
  499                 if self[i, j]:
  500                     return False
  501         return True
  502 
  503     def isLowerTriangularMatrix(self):
  504         """
  505         Check whether self is lower triangular matrix or not.
  506         """
  507         for i in range(1, self.row + 1):
  508             for j in range(i + 1, self.column + 1):
  509                 if self[i, j]:
  510                     return False
  511         return True
  512 
  513     def isDiagonalMatrix(self):
  514         """
  515         Check whether self is diagonal matrix or not.
  516         """
  517         return self.isUpperTriangularMatrix() and self.isLowerTriangularMatrix()
  518 
  519     def isScalarMatrix(self):
  520         """
  521         Check whether self is scalar matrix or not.
  522         Scalar matrix is matrix which is unit matrix times some scalar.
  523         """
  524         if not(self.isDiagonalMatrix()):
  525             return False
  526         chk = self[1, 1]
  527         for i in range(2, self.row + 1):
  528             if self[i, i] != chk:
  529                 return False
  530         return True
  531 
  532     def isSymmetricMatrix(self):
  533         """
  534         Check whether self is symmetric matrix or not.
  535         """
  536         for i in range(1, self.row + 1):
  537             for j in range(i + 1, self.column + 1):
  538                 if self[i, j] != self[j, i]:
  539                     return False
  540         return True
  541 
  542 
  543 class RingMatrix(Matrix):
  544     """
  545     RingMatrix is a class for matrices whose elements are in ring.
  546     """
  547 
  548     def __init__(self, row, column, compo=0, coeff_ring=0):
  549         """
  550         RingMatrix(row, column [,components, coeff_ring])
  551         """
  552         self._initialize(row, column, compo, coeff_ring)
  553         self._selectMatrix()
  554 
  555     def _selectMatrix(self):
  556         """
  557         Select Matrix class.
  558         """
  559         if self.row == self.column:
  560             self.__class__ = RingSquareMatrix
  561         else:
  562             self.__class__ = RingMatrix
  563 
  564     def __add__(self, other):
  565         """
  566         Return matrix addition.
  567         """
  568         if (self.row != other.row) or (self.column != other.column):
  569             raise MatrixSizeError()
  570         sums = []
  571         for i in range(1, self.row + 1):
  572             for j in range(1, self.column + 1):
  573                 sums.append(self[i, j] + other[i, j])
  574         return createMatrix(self.row, self.column, sums,
  575                   self.coeff_ring.getCommonSuperring(other.coeff_ring))
  576 
  577     def __sub__(self, other):
  578         """
  579         Return matrix subtraction.
  580         """
  581         if (self.row != other.row) or (self.column != other.column):
  582             raise MatrixSizeError()
  583         diff = []
  584         for i in range(1, self.row + 1):
  585             for j in range(1, self.column + 1):
  586                 diff.append(self[i, j] - other[i, j])
  587         return createMatrix(self.row, self.column, diff,
  588                   self.coeff_ring.getCommonSuperring(other.coeff_ring))
  589 
  590     def __mul__(self, other):
  591         """
  592         multiplication with a Matrix or a scalar
  593         """
  594         if isinstance(other, Matrix):
  595             if self.column != other.row:
  596                 raise MatrixSizeError()
  597             product = []
  598             for i in range(1, self.row + 1):
  599                 for j in range(1, other.column + 1):
  600                     part_product = self[i, 1] * other[1, j]
  601                     for k in range(2, self.column + 1):
  602                         part_product = part_product + self[i, k] * other[k, j]
  603                     product.append(part_product)
  604             return createMatrix(self.row, other.column, product,
  605                       self.coeff_ring.getCommonSuperring(other.coeff_ring))
  606         elif isinstance(other, vector.Vector):
  607             if self.column != len(other):
  608                 raise vector.VectorSizeError()
  609             product = []
  610             for i in range(1, self.row + 1):
  611                 part_product = self[i, 1] * other[1]
  612                 for j in range(2, self.column + 1):
  613                     part_product = part_product + self[i, j] * other[j]
  614                 product.append(part_product)
  615             return vector.Vector(product)
  616         else: #scalar mul
  617             try:
  618                 rings = self.coeff_ring.getCommonSuperring(ring.getRing(other))
  619             except:
  620                 return NotImplemented
  621             product = []
  622             for i in range(1, self.row + 1):
  623                 for j in range(1, self.column + 1):
  624                     product.append(self[i, j] * other)
  625             return createMatrix(self.row, self.column, product, rings)
  626 
  627     def __rmul__(self, other):
  628         if isinstance(other, Matrix):
  629             if other.column != self.row:
  630                 raise MatrixSizeError()
  631             product = []
  632             for i in range(1, other.row + 1):
  633                 for j in range(1, self.column + 1):
  634                     part_product = other[i, 1] * self[1, j]
  635                     for k in range(2, self.row + 1):
  636                         part_product = part_product + other[i, k] * self[k, j]
  637                     product.append(part_product)
  638             return createMatrix(other.row, self.column, product,
  639                      self.coeff_ring.getCommonSuperring(other.coeff_ring))
  640         elif isinstance(other, vector.Vector):
  641             if self.row != len(other):
  642                 raise vector.VectorSizeError()
  643             product = []
  644             for j in range(1, self.column + 1):
  645                 part_product = other[1] * self[1, j]
  646                 for i in range(2, self.row + 1):
  647                     part_product = part_product + other[i] * self[i, j]
  648                 product.append(part_product)
  649             return vector.Vector(product)
  650         else:
  651             try:
  652                 rings = self.coeff_ring.getCommonSuperring(ring.getRing(other))
  653             except:
  654                 return NotImplemented
  655             product = []
  656             for i in range(1, self.row + 1):
  657                 for j in range(1, self.column + 1):
  658                     product.append(self[i, j] * other)
  659             return createMatrix(self.row, self.column, product, rings)
  660 
  661     def __mod__(self, other):
  662         """
  663         return self modulo other.
  664         """
  665         if not bool(other):
  666             raise ZeroDivisionError()
  667         mod = []
  668         for i in range(1, self.row + 1):
  669             for j in range(1, self.column + 1):
  670                 mod.append(self[i, j] % other)
  671         return createMatrix(self.row, self.column, mod, self.coeff_ring)
  672 
  673     def __pos__(self):
  674         """
  675         return copy of self.
  676         """
  677         return self.copy()
  678 
  679     def __neg__(self):
  680         """
  681         return -self.
  682         """
  683         one = self.coeff_ring.one
  684         try:
  685             minus_one = -one
  686         except:
  687             minus_one = self.coeff_ring.zero - one
  688         return self.map(lambda ele: minus_one * ele)
  689 
  690     def getCoefficientRing(self):
  691         """
  692         Set and return coefficient ring.
  693         """
  694         if not hasattr(self, "_coeff_ring"):
  695             scalars = None
  696             for i in range(1, self.row + 1):
  697                 for j in range(1, self.column + 1):
  698                     cring = ring.getRing(self[i, j])
  699                     if scalars is None or \
  700                     scalars != cring and scalars.issubring(cring):
  701                         scalars = cring
  702                     elif not scalars.issuperring(cring):
  703                         scalars = scalars.getCommonSuperring(cring)
  704             if scalars != self.coeff_ring or scalars.issubring(self.coeff_ring):
  705                 scalars = self.coeff_ring
  706             elif not scalars.issuperring(self.coeff_ring):
  707                 scalars = scalars.getCommonSuperring(self.coeff_ring)
  708             self._coeff_ring = self.coeff_ring = scalars
  709         return self._coeff_ring
  710 
  711     def toFieldMatrix(self):
  712         """RingMatrix -> FieldMatrix"""
  713         self.__class__ = FieldMatrix
  714         self.coeff_ring = self.coeff_ring.getQuotientField()
  715 
  716     def toSubspace(self, isbasis=None):
  717         """RingMatrix -> Subspace"""
  718         self.toFieldMatrix()
  719         self.toSubspace(isbasis)
  720 
  721     def hermiteNormalForm(self, non_zero=False):
  722         # Algorithm 2.4.5 in CCANT (fixed for parameter l)
  723         """
  724         Find the Hermite Normal Form for integer matrix.
  725         If non_zero is True, return non-zero columns for self's HNF
  726         """
  727         A = self.copy()
  728         rings = self.coeff_ring
  729         # step 1 [Initialize]
  730         j0 = k = self.column
  731         for i in range(1, self.row + 1)[::-1]:
  732             while 1:
  733                 # step 2 [Check zero]
  734                 for j in range(1, j0)[::-1]:
  735                     if bool(A[i, j]):
  736                         break
  737                 else: # j==1
  738                     break
  739                 # step 3 [Euclidean step]
  740                 u, v, d = rings.extgcd(A[i, k], A[i, j])
  741                 A_ik = ring.exact_division(A[i, k], d)
  742                 A_ij = ring.exact_division(A[i, j], d)
  743                 B = u * A[k] + v * A[j]
  744                 A[j] = A_ik * A[j] - A_ij * A[k]
  745                 A[k] = B
  746             # step4 [Final reductions]
  747             b = A[i, k]
  748             if b < 0:
  749                 A[k] = -A[k]
  750                 b = -b
  751             if not bool(b):
  752                 k += 1
  753             else:
  754                 for j in range(k + 1, self.column + 1):
  755                     q = A[i, j] // b
  756                     A[j] = A[j] - q * A[k]
  757             # step 5 [Finished?]
  758             k -= 1
  759             j0 = k
  760             if k == 0:
  761                 break
  762             # go to step 2
  763         if non_zero:
  764             if k == self.column: # zero module
  765                 return None
  766             W = createMatrix(
  767                 self.row, self.column - k,
  768                 [A[j + k] for j in range(1, self.column - k + 1)])
  769             return W
  770         else:
  771             return A
  772 
  773     HNF = hermiteNormalForm
  774 
  775     def _SimplifyHNF(self):
  776         """
  777         This method is a common process used in 
  778         extHNF() and kernelAsModule()
  779         """
  780         A = self.copy()
  781         U = unitMatrix(A.column, A.coeff_ring)
  782         rings = self.coeff_ring
  783         # step 1 [Initialize]
  784         j0 = k = self.column
  785         for i in range(1, self.row + 1)[::-1]:
  786             while 1:
  787                 # step 2 [Check zero]
  788                 for j in range(1, j0)[::-1]:
  789                     if bool(A[i, j]):
  790                         break
  791                 else: # j==1
  792                     break
  793                 # step 3 [Euclidean step]
  794                 u, v, d = rings.extgcd(A[i, k], A[i, j])
  795                 A_ik = ring.exact_division(A[i, k], d)
  796                 A_ij = ring.exact_division(A[i, j], d)
  797                 B = u * A[k] + v * A[j]
  798                 A[j] = A_ik * A[j] - A_ij * A[k]
  799                 A[k] = B
  800                 B = u * U[k] + v * U[j]
  801                 U[j] = A_ik * U[j] - A_ij * U[k]
  802                 U[k] = B
  803             # step4 [Final reductions]
  804             b = A[i, k]
  805             if b < 0:
  806                 A[k] = -A[k]
  807                 U[k] = -U[k]
  808                 b = -b
  809             if not bool(b):
  810                 k += 1
  811             else:
  812                 for j in range(k + 1, self.column + 1):
  813                     q = A[i, j] // b
  814                     A[j] = A[j] - q * A[k]
  815                     U[j] = U[j] - q * U[k]
  816             # step 5 [Finished?]
  817             k -= 1
  818             j0 = k
  819             if k == 0:
  820                 break
  821             # go to step 2
  822         return (U, A, k)
  823 
  824     def exthermiteNormalForm(self, non_zero=False):
  825         # Modified Algorithm 2.4.5 in CCANT
  826         """
  827         Find the Hermite Normal Form M for integer matrix.
  828         Computing U which satisfied M=self*U.
  829         Return matrices tuple,(U, M).
  830         """
  831         U, A, k = self._SimplifyHNF()
  832         if non_zero:
  833             if k == self.column: # zero module
  834                 return None
  835             new_A = createMatrix(
  836                 self.row, self.column - k,
  837                 [A[j + k] for j in range(1, self.column - k + 1)])
  838             new_U = createMatrix(
  839                 self.column, self.column - k,
  840                 [U[j + k] for j in range(1, self.column - k + 1)])
  841             return (new_U, new_A)
  842         else:
  843             return (U, A)
  844 
  845     extHNF = exthermiteNormalForm
  846 
  847     def kernelAsModule(self):
  848         """
  849         Compute kernel as Z-module.
  850         """
  851         U, A, k = self._SimplifyHNF()
  852         if k == 0:
  853             return None
  854         else:
  855             ker = createMatrix(
  856                 self.column, k, [U[j] for j in range(1, k + 1)])
  857             return ker
  858 
  859 
  860 class RingSquareMatrix(SquareMatrix, RingMatrix, ring.RingElement):
  861     """
  862     RingSquareMatrix is a class for square matrices whose elements are in ring.
  863     """
  864 
  865     def __init__(self, row, column=0, compo=0, coeff_ring=0):
  866         """
  867         RingSquareMatrix(row [, column ,components, coeff_ring])
  868         RingSquareMatrix must be row == column .
  869         """
  870         self._initialize(row, column, compo, coeff_ring)
  871 
  872     def __pow__(self, other):
  873         """
  874         powering self to integer.
  875         """
  876         n = +other
  877         if not isinstance(n, (int, long)):
  878             raise TypeError("index must be an integer")
  879         power = unitMatrix(self.row, self.coeff_ring)
  880         # check n
  881         if n == 0:
  882             return power
  883         if n > 0:
  884             z = self.copy()
  885         else:
  886             if hasattr(self, "inverse"):
  887                 n = abs(n)
  888                 z = self.inverse()
  889             else:
  890                 raise NoInverse()
  891         while 1:
  892             if n & 1:
  893                 power = power * z
  894             n //= 2
  895             if n == 0:
  896                 return power
  897             z = z * z
  898 
  899     def toFieldMatrix(self):
  900         """RingSquareMatrix -> FieldSquareMatrix"""
  901         self.__class__ = FieldSquareMatrix
  902         self.coeff_ring = self.coeff_ring.getQuotientField()
  903 
  904     def getRing(self):
  905         """
  906         Return matrix ring of self.
  907         """
  908         return MatrixRing.getInstance(self.row, self.getCoefficientRing())
  909 
  910     def isOrthogonalMatrix(self):
  911         """
  912         Check whether self is orthogonal matrix or not.
  913         Orthogonal matrix satisfies M*M^T equals unit matrix.
  914         """
  915         return self * self.transpose() == unitMatrix(self.row, self.coeff_ring)
  916 
  917     def isAlternatingMatrix(self):
  918         """
  919         Check whether self is alternating matrix or not.
  920         Alternating (skew symmetric, or antisymmetric) matrix satisfies M=-M^T.
  921         """
  922         for i in range(1, self.row + 1):
  923             for j in range(i, self.column + 1):
  924                 if self[i, j] != -self[j, i]:
  925                     return False
  926         return True
  927 
  928     isAntisymmetricMatrix = isAlternatingMatrix
  929     isSkewsymmetricMatrix = isAlternatingMatrix
  930 
  931     def isSingular(self):
  932         """
  933         Check determinant == 0 or not.
  934         """
  935         return not bool(self.determinant())
  936 
  937     def trace(self):
  938         """
  939         Return trace of self.
  940         """
  941         trace = self.coeff_ring.zero
  942         for i in range(1, self.row + 1):
  943             trace = trace + self[i, i]
  944         return trace
  945 
  946     def determinant(self): # Algorithm 2.2.6 of Cohen's book
  947         """
  948         Return determinant of self.
  949         """
  950         M = self.copy()
  951         n = self.row
  952         c = self.coeff_ring.one
  953         sign = True
  954         for k in range(1, n):
  955             p = M[k, k]
  956             if not bool(p): # p==0
  957                 i = k + 1
  958                 while not bool(M[i, k]):
  959                     if i == n:
  960                         return self.coeff_ring.zero
  961                     else:
  962                         i += 1
  963                 for j in range(k, n + 1):
  964                     tmp = M[i, j]
  965                     M[i, j] = M[k, j]
  966                     M[k, j] = tmp
  967                 sign = not(sign)
  968                 p = M[k, k]
  969             for i in range(k + 1, n + 1):
  970                 for j in range(k + 1, n + 1):
  971                     t = p * M[i, j] - M[i, k] * M[k, j]
  972                     M[i, j] = ring.exact_division(t, c)
  973             c = p
  974         if sign:
  975             return M[n, n]
  976         else:
  977             return -M[n, n]
  978 
  979     def cofactor(self, i, j):
  980         """
  981         Return (i, j)-cofactor of self.
  982         """
  983         cofactor = (self.subMatrix(i, j)).determinant()
  984         if (i + j) & 1:
  985             cofactor = -cofactor
  986         return cofactor
  987 
  988     def commutator(self, other):
  989         """
  990         Return commutator defined as follows:
  991         [self, other] = self * other - other * self .
  992         """
  993         return self * other - other * self
  994 
  995     def characteristicMatrix(self):
  996         """
  997         Return the characteristic matrix (i.e. xI-A) of self.
  998         """
  999         import nzmath.poly.uniutil as uniutil
 1000         x = uniutil.polynomial({1:1}, self.coeff_ring)
 1001         return x * unitMatrix(self.row, x.getRing()) - self
 1002 
 1003     def _characteristicPolyList(self): # Algorithm 2.2.7 of Cohen's book
 1004         """
 1005         for characteristicPolynomial, adjugateMatrix
 1006         
 1007         Assume self.row >= 2.
 1008         """
 1009         unit = unitMatrix(self.row, self.coeff_ring)
 1010         coeff = [self.coeff_ring.one, -self.trace()]
 1011         C = self + coeff[-1] * unit
 1012         i = 2
 1013         while i < self.row:
 1014             C = self * C
 1015             coeff.append(ring.exact_division(-C.trace(), i))
 1016             C = C + coeff[-1] * unit
 1017             i += 1
 1018         coeff.append(ring.exact_division(-(self * C).trace(), i))
 1019         coeff.reverse()
 1020         return coeff, C
 1021 
 1022     def characteristicPolynomial(self):
 1023         """
 1024         characteristicPolynomial() -> Polynomial
 1025         """
 1026         import nzmath.poly.uniutil
 1027         genPoly = nzmath.poly.uniutil.polynomial
 1028         if self.row == 1:
 1029             rings = self.coeff_ring
 1030             return genPoly({0:-self.trace(), 1:rings.one}, rings)
 1031         coeff = self._characteristicPolyList()[0]
 1032         return genPoly(dict(enumerate(coeff)), self.coeff_ring)
 1033 
 1034     def adjugateMatrix(self):
 1035         """
 1036         Return adjugate(classical adjoint) matrix.
 1037         """
 1038         if self.row == 1:
 1039             return unitMatrix(self.row, self.coeff_ring)
 1040         C = self._characteristicPolyList()[1]
 1041         if self.row & 1:
 1042             return C
 1043         else:
 1044             return -C
 1045 
 1046     def cofactorMatrix(self):
 1047         """
 1048         Return cofactor matrix.
 1049         """
 1050         return self.adjugateMatrix().transpose()
 1051 
 1052     cofactors = cofactorMatrix
 1053 
 1054     def smithNormalForm(self):# Algorithm 2.4.14 of Cohen's book
 1055         """
 1056         Find the Smith Normal Form for square non-singular integral matrix.
 1057         Return the list of diagonal elements.
 1058         """
 1059         M = self.copy()
 1060         n = M.row
 1061         R = M.determinant()
 1062         rings = self.coeff_ring
 1063         if not bool(R):
 1064             raise ValueError("Don't input singular matrix")
 1065         if R < 0:
 1066             R = -R
 1067         lst = []
 1068         while n != 1:
 1069             j = n
 1070             c = 0
 1071             while j != 1:
 1072                 j -= 1
 1073                 if M[n, j]:
 1074                     u, v, d = rings.extgcd(M[n, j], M[n, n])
 1075                     B = v * M.getColumn(n) + u * M.getColumn(j)
 1076                     M_nn = ring.exact_division(M[n, n], d)
 1077                     M_nj = ring.exact_division(M[n, j], d)
 1078                     M.setColumn(j, ((M_nn * M.getColumn(j)
 1079                                      - M_nj * M.getColumn(n)) % R))
 1080                     M.setColumn(n, (B % R))
 1081             j = n
 1082             while j != 1:
 1083                 j -= 1
 1084                 if M[j, n]:
 1085                     u, v, d = rings.extgcd(M[j, n], M[n, n])
 1086                     B = v * M.getRow(n) + u * M.getRow(j)
 1087                     M_nn = ring.exact_division(M[n, n], d)
 1088                     M_jn = ring.exact_division(M[j, n], d)
 1089                     M.setRow(j, ((M_nn * M.getRow(j)
 1090                                   - M_jn * M.getRow(n)) % R))
 1091                     M.setRow(n, (B % R))
 1092                     c += 1
 1093             if c <= 0:
 1094                 b = M[n, n]
 1095                 flag = False
 1096                 if not bool(b):
 1097                     b = R
 1098                 for k in range(1, n):
 1099                     for l in range(1, n):
 1100                         if (M[k, l] % b):
 1101                             M.setRow(n, M.getRow(n) + M.getRow(k))
 1102                             flag = True
 1103                 if not flag:
 1104                     dd = rings.gcd(M[n, n], R)
 1105                     lst.append(dd)
 1106                     R = ring.exact_division(R, dd)
 1107                     n -= 1
 1108         dd = rings.gcd(M[1, 1], R)
 1109         lst.append(dd)
 1110         lst.reverse()
 1111         return lst
 1112 
 1113     SNF = smithNormalForm
 1114     elementary_divisor = smithNormalForm
 1115 
 1116     def extsmithNormalForm(self):
 1117         """
 1118         Find the Smith Normal Form M for square matrix,
 1119         Computing U,V which satisfied M=U*self*V.
 1120         Return matrices tuple,(U,V,M).
 1121         """
 1122         M = self.copy()
 1123         n = M.row
 1124         U = unitMatrix(M.row, M.coeff_ring)
 1125         V = unitMatrix(M.row, M.coeff_ring)
 1126         rings = self.coeff_ring
 1127         while n != 1:
 1128             j = n
 1129             c = 0
 1130             while j != 1:
 1131                 j -= 1
 1132                 if M[n, j]:
 1133                     u, v, d = rings.extgcd(M[n, j], M[n, n])
 1134                     M_nn = ring.exact_division(M[n, n], d)
 1135                     M_nj = ring.exact_division(M[n, j], d)
 1136                     B = v * M.getColumn(n) + u * M.getColumn(j)
 1137                     M.setColumn(j, (M_nn * M.getColumn(j)
 1138                                     - M_nj * M.getColumn(n)))
 1139                     M.setColumn(n, B)
 1140                     B = v * V.getColumn(n) + u * V.getColumn(j)
 1141                     V.setColumn(j, (M_nn * V.getColumn(j)
 1142                                     - M_nj * V.getColumn(n)))
 1143                     V.setColumn(n, B)
 1144             j = n
 1145             while j != 1:
 1146                 j -= 1
 1147                 if M[j, n]:
 1148                     u, v, d = rings.extgcd(M[j, n], M[n, n])
 1149                     M_nn = ring.exact_division(M[n, n], d)
 1150                     M_jn = ring.exact_division(M[j, n], d)
 1151                     B = v * M.getRow(n) + u * M.getRow(j)
 1152                     M.setRow(j, (M_nn * M.getRow(j) - M_jn * M.getRow(n)))
 1153                     M.setRow(n, B)
 1154                     B = v * U.getRow(n) + u * U.getRow(j)
 1155                     U.setRow(j, (M_nn * U.getRow(j) - M_jn * U.getRow(n)))
 1156                     U.setRow(n, B)
 1157                     c += 1
 1158             if c <= 0:
 1159                 b = M[n, n]
 1160                 flag = False
 1161                 for k in range(1, n):
 1162                     for l in range(1, n):
 1163                         if (M[k, l] % b):
 1164                             M.setRow(n, M.getRow(n) + M.getRow(k))
 1165                             U.setRow(n, U.getRow(n) + U.getRow(k))
 1166                             flag = True
 1167                 if not flag:
 1168                     n -= 1
 1169         for j in range(1, M.column + 1):
 1170             if M[j, j] < 0:
 1171                 V[j] = -V[j]
 1172                 M[j, j] = -M[j, j]
 1173         return (U, V, M)
 1174 
 1175     extSNF = extsmithNormalForm
 1176 
 1177 
 1178 class FieldMatrix(RingMatrix):
 1179     """
 1180     FieldMatrix is a class for matrices whose elements are in field.
 1181     """
 1182 
 1183     def __init__(self, row, column, compo=0, coeff_ring=0):
 1184         """
 1185         FieldMatrix(row, column [,components, coeff_ring])
 1186         """
 1187         self._initialize(row, column, compo, coeff_ring)
 1188         self._selectMatrix()
 1189 
 1190     def _initialize(self, row, column, compo=0, coeff_ring=0):
 1191         """initialize matrix"""
 1192         RingMatrix._initialize(self, row, column, compo, coeff_ring)
 1193         if not self.coeff_ring.isfield():
 1194             self.coeff_ring = self.coeff_ring.getQuotientField()
 1195 
 1196     def _selectMatrix(self):
 1197         """
 1198         Select Matrix class.
 1199         """
 1200         if self.__class__ != Subspace:
 1201             if self.row == self.column:
 1202                 self.__class__ = FieldSquareMatrix
 1203             else:
 1204                 self.__class__ = FieldMatrix
 1205 
 1206     def __truediv__(self, other):
 1207         """
 1208         division by a scalar.
 1209         """
 1210         return ring.inverse(other) * self
 1211 
 1212     __div__ = __truediv__ # backward compatibility?
 1213 
 1214     def toSubspace(self, isbasis=None):
 1215         """FieldMatrix -> Subspace"""
 1216         self.__class__ = Subspace
 1217         self.isbasis = isbasis
 1218 
 1219     def _cohensSimplify(self):
 1220         """
 1221         _cohensSimplify is a common process used in image() and kernel()
 1222 
 1223         Return a tuple of modified matrix M, image data c and kernel data d.
 1224         """
 1225         M = self.copy()
 1226         c = [0] * (M.row + 1)
 1227         d = [-1] * (M.column + 1)
 1228         for k in range(1, M.column + 1):
 1229             for j in range(1, M.row + 1):
 1230                 if not c[j] and M[j, k]:
 1231                     break
 1232             else:           # not found j such that m(j, k)!=0 and c[j]==0
 1233                 d[k] = 0
 1234                 continue
 1235             top = -ring.inverse(M[j, k])
 1236             M[j, k] = -self.coeff_ring.one
 1237             for s in range(k + 1, M.column + 1):
 1238                 M[j, s] = top * M[j, s]
 1239             for i in range(1, M.row + 1):
 1240                 if i == j:
 1241                     continue
 1242                 top = M[i, k]
 1243                 M[i, k] = self.coeff_ring.zero
 1244                 for s in range(k + 1, M.column + 1):
 1245                     M[i, s] = M[i, s] + top * M[j, s]
 1246             c[j] = k
 1247             d[k] = j
 1248         return (M, c, d)
 1249 
 1250     def kernel(self):       # Algorithm 2.3.1 of Cohen's book
 1251         """
 1252         Return a Matrix whose column vectors are one basis of self's kernel,
 1253         or return None if self's kernel is 0.
 1254         """
 1255         tmp = self._cohensSimplify()
 1256         M, d = tmp[0], tmp[2]
 1257         basis = []
 1258         for k in range(1, M.column + 1):
 1259             if d[k]:
 1260                 continue
 1261             vect = []
 1262             for i in range(1, M.column + 1):
 1263                 if d[i] > 0:
 1264                     vect.append(M[d[i], k])
 1265                 elif i == k:
 1266                     vect.append(self.coeff_ring.one)
 1267                 else:
 1268                     vect.append(self.coeff_ring.zero)
 1269             basis.append(vect)
 1270         dimension = len(basis)
 1271         if dimension == 0:
 1272             return None
 1273         output = zeroMatrix(self.column, dimension, self.coeff_ring)
 1274         for j in range(1, dimension + 1):
 1275             output.setColumn(j, basis[j - 1])
 1276         return output
 1277 
 1278     def image(self):        # Algorithm 2.3.2 of Cohen's book
 1279         """
 1280         Return a Matrix which column vectors are one basis of self's image,
 1281         or return None if self's image is 0.
 1282         """
 1283         tmp = self._cohensSimplify()
 1284         M, c = tmp[0], tmp[1]
 1285         basis = []
 1286         for j in range(1, M.row + 1):
 1287             if c[j]:
 1288                 basis.append(self[c[j]])
 1289         dimension = len(basis)
 1290         if dimension == 0:
 1291             return None
 1292         output = zeroMatrix(self.row, dimension, self.coeff_ring)
 1293         for j in range(1, dimension + 1):
 1294             output.setColumn(j, basis[j - 1])
 1295         output._selectMatrix()
 1296         return output
 1297 
 1298     def rank(self):
 1299         """
 1300         Return rank of self.
 1301         """
 1302         img = self.image()
 1303         if img:
 1304             return img.column
 1305         else:
 1306             return 0
 1307 
 1308     def inverseImage(self, V):    # modified Algorithm 2.3.5 of Cohen's book
 1309         """
 1310         inverseImage(V) -> X
 1311         
 1312         such that self * X == V
 1313         """
 1314         if isinstance(V, vector.Vector):
 1315             if self.row != len(V):
 1316                 raise vector.VectorSizeError()
 1317             B = createMatrix(len(V), 1, V.compo)
 1318         else:
 1319             if self.row != V.row:
 1320                 raise MatrixSizeError()
 1321             B = V.copy() # step 1
 1322         M = self.copy()
 1323         m = M.row
 1324         n = M.column
 1325         r = B.column
 1326         X = zeroMatrix(n, r, self.coeff_ring)
 1327         non_zero = []
 1328         i = 1
 1329         # step 2
 1330         for j in range(1, n + 1):
 1331             # step 3
 1332             for k in range(i, m + 1):
 1333                 if M[k, j]:
 1334                     break
 1335             else:
 1336                 continue
 1337             # step 4
 1338             if k > i:
 1339                 for l in range(j, n + 1):
 1340                     t = M[i, l]
 1341                     M[i, l] = M[k, l]
 1342                     M[k, l] = t
 1343                 B.swapRow(i, k)
 1344             # step 5
 1345             d = ring.inverse(M[i, j])
 1346             for k in range(i + 1, m + 1):
 1347                 ck = d * M[k, j]
 1348                 for l in range(j + 1, n + 1):
 1349                     M[k, l] = M[k, l] - ck * M[i, l]
 1350                 for l in range(r):
 1351                     B[k, l] = B[k, l] - ck * B[i, l]
 1352             non_zero.insert(0, j)
 1353             i += 1
 1354             if i > m:
 1355                 break
 1356         # step 6
 1357         i -= 1
 1358         zero = self.coeff_ring.zero
 1359         for j in non_zero:
 1360             d = ring.inverse(M[i, j])
 1361             for k in range(r):
 1362                 sums = zero
 1363                 for l in range(j + 1, n + 1):
 1364                     sums = sums + M[i, l] * X[l, k]
 1365                 X[j, k] = (B[i, k] - sums) * d
 1366             i -= 1
 1367         # step 7
 1368         i = len(non_zero) + 1
 1369         for j in range(1, r + 1):
 1370             for k in range(i, m + 1):
 1371                 if B[k, j]:
 1372                     raise NoInverseImage()
 1373         return X
 1374 
 1375     def solve(self, B):  # modified Algorithm 2.3.4 of Cohen's book
 1376         """
 1377         Return solution X for self * X = B (B is a vector).
 1378         This function returns tuple (V, M) below.
 1379           V: one solution as vector
 1380           M: kernel of self as list of basis vectors.
 1381         If you want only one solution, use 'inverseImage'.
 1382 
 1383         Warning: B should not be a matrix instead of a vector
 1384         """
 1385         M_1 = self.copy()
 1386         M_1.insertColumn(self.column + 1, B.compo)
 1387         V = M_1.kernel()
 1388         ker = []
 1389         flag = False
 1390         if not V:
 1391             raise NoInverseImage("no solution")
 1392         n = V.row
 1393         for j in range(1, V.column + 1):
 1394             if not bool(V[n, j]): # self's kernel
 1395                 ker.append(vector.Vector([V[i, j] for i in range(1, n)]))
 1396             elif not(flag):
 1397                 d = -ring.inverse(V[n, j])
 1398                 sol = vector.Vector([V[i, j] * d for i in range(1, n)])
 1399                 flag = True
 1400         if not(flag):
 1401             raise NoInverseImage("no solution")
 1402         return sol, ker
 1403 
 1404     def columnEchelonForm(self):  # Algorithm 2.3.11 of Cohen's book
 1405         """
 1406         Return a Matrix in column echelon form whose image is equal to 
 1407         the image of self.
 1408         """
 1409         M = self.copy()
 1410         k = M.column
 1411         for i in range(M.row, 0, -1):
 1412             for j in range(k, 0, -1):
 1413                 if M[i, j]:
 1414                     break
 1415             else:
 1416                 continue
 1417             d = ring.inverse(M[i, j])
 1418             for l in range(1, i + 1):
 1419                 t = d * M[l, j]
 1420                 M[l, j] = M[l, k]
 1421                 M[l, k] = t
 1422             for j in range(1, M.column + 1):
 1423                 if j == k:
 1424                     continue
 1425                 for l in range(1, i + 1):
 1426                     M[l, j] = M[l, j] - M[l, k] * M[i, j]
 1427             k -= 1
 1428         return M
 1429 
 1430 
 1431 class FieldSquareMatrix(RingSquareMatrix, FieldMatrix):
 1432     """
 1433     FieldSquareMatrix is a class for square matrices in field.
 1434     """
 1435 
 1436     def __init__(self, row, column=0, compo=0, coeff_ring=0):
 1437         """
 1438         FieldSquareMatrix(row [, column, components, coeff_ring])
 1439         FieldSquareMatrix must be row == column .
 1440         """
 1441         self._initialize(row, column, compo, coeff_ring)
 1442 
 1443     def triangulate(self):
 1444         """
 1445         Return triangulated matrix of self.
 1446         """
 1447         triangle = self.copy()
 1448         flag = False # for calculation of determinant
 1449         for i in range(1, triangle.row + 1):
 1450             if not triangle[i, i]:
 1451                 for k in range(i + 1, triangle.row + 1):
 1452                     if triangle[k, i]:
 1453                         triangle.swapRow(i + 1, k + 1)
 1454                         flag = not(flag)
 1455                         break # break the second loop
 1456                 else:
 1457                     # the below components are all 0. Back to the first loop
 1458                     continue
 1459             for k in range(i + 1, triangle.row + 1):
 1460                 inv_i_i = ring.inverse(triangle[i, i])
 1461                 ratio = triangle[k, i] * inv_i_i
 1462                 for l in range(i, triangle.column + 1):
 1463                     triangle[k, l] = triangle[k, l] - triangle[i, l] * ratio
 1464         if flag:
 1465             for j in range(triangle.row, triangle.column + 1):
 1466                 triangle[triangle.row, j] = -triangle[triangle.row, j]
 1467         return triangle
 1468 
 1469     def determinant(self):
 1470         """
 1471         Return determinant of self.
 1472         """
 1473         triangle = self.triangulate()
 1474         det = self.coeff_ring.one
 1475         for i in range(1, self.row + 1):
 1476             det = det * triangle[i, i]
 1477         return det
 1478 
 1479     def inverse(self, V=1): # modified Algorithm 2.2.2, 2.3.5 of Cohen's book
 1480         """
 1481         Return inverse matrix of self or self.inverse() * V.
 1482         If inverse does not exist, raise NoInverse error.
 1483         """
 1484         if isinstance(V, vector.Vector):
 1485             if self.row != len(V):
 1486                 raise vector.VectorSizeError()
 1487             B = createMatrix(len(V), 1, V.compo)
 1488         elif isinstance(V, Matrix):
 1489             if self.row != V.row:
 1490                 raise MatrixSizeError()
 1491             B = V.copy() # step 1
 1492         else: # V==1
 1493             B = unitMatrix(self.row, self.coeff_ring)
 1494         M = self.copy()
 1495         n = M.row
 1496         r = B.column
 1497         X = zeroMatrix(n, r, self.coeff_ring)
 1498         # step 2
 1499         for j in range(1, n + 1):
 1500             # step 3
 1501             for i in range(j, n + 1):
 1502                 if M[i, j]:
 1503                     break
 1504             else:
 1505                 raise NoInverse()
 1506             # step 4
 1507             if i > j:
 1508                 for l in range(j, n + 1):
 1509                     t = M[i, l]
 1510                     M[i, l] = M[j, l]
 1511                     M[j, l] = t
 1512                 B.swapRow(i, j)
 1513             # step 5
 1514             d = ring.inverse(M[j, j])
 1515             for k in range(j + 1, n + 1):
 1516                 ck = d * M[k, j]
 1517                 for l in range(j + 1, n + 1):
 1518                     M[k, l] = M[k, l] - ck * M[j, l]
 1519                 for l in range(1, r + 1):
 1520                     B[k, l] = B[k, l] - ck * B[j, l]
 1521         # step 6
 1522         for i in range(n, 0, -1):
 1523             d = ring.inverse(M[i, i])
 1524             for k in range(1, r + 1):
 1525                 sums = self.coeff_ring.zero
 1526                 for j in range(i + 1, n + 1):
 1527                     sums = sums + M[i, j] * X[j, k]
 1528                 X[i, k] = (B[i, k] - sums) * d
 1529         if r != 1:
 1530             return X
 1531         else:
 1532             return X[1]
 1533 
 1534     def hessenbergForm(self):      # Algorithm 2.2.9 of Cohen's book
 1535         """Return a Matrix in Hessenberg Form."""
 1536         n = self.row
 1537         zero = self.coeff_ring.zero
 1538         # step 1
 1539         H = self.copy()
 1540         for m in range(2, H.row):
 1541             # step 2
 1542             for i in range(m + 1, n + 1):
 1543                 if H[i, m - 1]:
 1544                     break
 1545             else:
 1546                 continue
 1547             tinv = ring.inverse(H[i, m - 1])
 1548             if i > m:
 1549                 for j in range(m - 1, n + 1):
 1550                     tmp = H[i, j]
 1551                     H[i, j] = H[m, j]
 1552                     H[m, j] = tmp
 1553                 H.swapColumn(i, m)
 1554             # step 3
 1555             for i in range(m + 1, n + 1):
 1556                 if H[i, m - 1]:
 1557                     u = H[i, m - 1] * tinv
 1558                     for j in range(m, n + 1):
 1559                         H[i, j] = H[i, j] - u * H[m, j]
 1560                     H[i, m - 1] = zero
 1561                     H.setColumn(m, H[m] + u * H[i])
 1562         return H
 1563 
 1564     def LUDecomposition(self):
 1565         """
 1566         LUDecomposition() -> (L, U)
 1567         
 1568         L and U are matrices such that
 1569             self == L * U
 1570             L : lower triangular matrix
 1571             U : upper triangular matrix
 1572         """
 1573 
 1574         n = self.row
 1575         L = unitMatrix(n, self.coeff_ring)
 1576         U = self.copy()
 1577         # initialize L and U
 1578         for i in range(1, n + 1):
 1579             for j in range(i + 1, n + 1):
 1580                 L[j, i] = U[j, i] * ring.inverse(U[i, i])
 1581                 for k in range(i, n + 1):
 1582                     U[j, k] = U[j, k] - U[i, k] * L[j, i]
 1583         return (L, U)
 1584 
 1585 
 1586 class MatrixRing (ring.Ring):
 1587     """
 1588     MatrixRing is a class for matrix rings.
 1589     """
 1590 
 1591     _instances = {}
 1592 
 1593     def __init__(self, size, scalars):
 1594         """
 1595         MatrixRing(size, scalars)
 1596         
 1597         size: size of matrices (positive integer)
 1598         scalars: ring of scalars
 1599         """
 1600         ring.Ring.__init__(self)
 1601         self.size = size
 1602         self.scalars = scalars
 1603 
 1604     def __eq__(self, other):
 1605         """
 1606         self == other
 1607         """
 1608         return (self.__class__ == other.__class__ and
 1609                 self.size == other.size and
 1610                 self.scalars == other.scalars)
 1611 
 1612     def __hash__(self):
 1613         return self.scalars ** self.size
 1614 
 1615     def __repr__(self):
 1616         return "MatrixRing(%d, %s)" % (self.size, self.scalars)
 1617 
 1618     def __str__(self):
 1619         return "M_%d(%s)" % (self.size, str(self.scalars))
 1620 
 1621     @classmethod
 1622     def getInstance(cls, size, scalars):
 1623         """
 1624         Return the cached instance of the specified matrix ring.  If
 1625         the specified ring is not cached, it is created, cached and
 1626         returned.
 1627         
 1628         The method is a class method.
 1629         """
 1630         if (size, scalars) not in cls._instances:
 1631             anInstance = MatrixRing(size, scalars)
 1632             cls._instances[size, scalars] = anInstance
 1633         return cls._instances[size, scalars]
 1634 
 1635     def unitMatrix(self):
 1636         """
 1637         Return the unit matrix.
 1638         """
 1639         return self.one.copy()
 1640 
 1641     def _getOne(self):
 1642         """
 1643         getter for one (unit matrix)
 1644         """
 1645         if self._one is None:
 1646             self._one = unitMatrix(self.size, self.scalars)
 1647         return self._one
 1648 
 1649     one = property(_getOne, None, None, "multiplicative unit")
 1650 
 1651     def zeroMatrix(self):
 1652         """
 1653         Return the zero matrix.
 1654         """
 1655         return self.zero.copy()
 1656 
 1657     def _getZero(self):
 1658         """
 1659         Return zero matrix.
 1660         """
 1661         if self._zero is None:
 1662             self._zero = zeroMatrix(self.size, self.scalars)
 1663         return self._zero
 1664 
 1665     zero = property(_getZero, None, None, "additive unit")
 1666 
 1667     def createElement(self, compo):
 1668         """
 1669         Return a newly created matrix from 'compo'.
 1670 
 1671         'compo' must be a list of n*n components in the scalar ring,
 1672         where n = self.size.
 1673         """
 1674         return createMatrix(self.size, compo, self.scalars)
 1675 
 1676     def getCharacteristic(self):
 1677         """
 1678         Return the characteristic of the ring.
 1679         """
 1680         return self.scalars.getCharacteristic()
 1681 
 1682     def issubring(self, other):
 1683         """
 1684         Report whether another ring contains the ring as a subring.
 1685         """
 1686         if other is self:
 1687             return True
 1688         if not isinstance(other, MatrixRing):
 1689             return False
 1690         return self.size == other.size and self.scalars.issubring(other.scalars)
 1691 
 1692     def issuperring(self, other):
 1693         """
 1694         Report whether the ring is a superring of another ring.
 1695         """
 1696         if other is self:
 1697             return True
 1698         if not isinstance(other, MatrixRing):
 1699             return False
 1700         return self.size == other.size and \
 1701         self.scalars.issuperring(other.scalars)
 1702 
 1703     def getCommonSuperring(self, other):
 1704         """
 1705         Return common super ring of self and another ring.
 1706         """
 1707         if not isinstance(other, MatrixRing) or self.size != other.size:
 1708             raise TypeError("no common super ring")
 1709         return MatrixRing.getInstance(self.size, 
 1710         self.scalars.getCommonSuperring(other.scalars))
 1711 
 1712 
 1713 class Subspace(FieldMatrix):
 1714     """
 1715     Subspace is a class for subspaces.
 1716     """
 1717 
 1718     def __init__(self, row, column, compo=0, coeff_ring=0, isbasis=None):
 1719         """
 1720         Subspace(row, column [,components, coeff_ring, isbasis])
 1721         """
 1722         if isinstance(compo, bool):
 1723             isbasis = compo
 1724             compo = 0
 1725         elif isinstance(coeff_ring, bool):
 1726             isbasis = coeff_ring
 1727             coeff_ring = 0
 1728         self._initialize(row, column, compo, coeff_ring)
 1729         self.isbasis = isbasis
 1730 
 1731     @classmethod
 1732     def fromMatrix(cls, mat, isbasis=None):
 1733         """
 1734         A constructor class method, which creates Subspace from a
 1735         Matrix instance.
 1736         """
 1737         compo = []
 1738         for row in mat.compo:
 1739             compo += row
 1740         return cls(mat.row, mat.column, compo, mat.coeff_ring, isbasis)
 1741 
 1742     def toFieldMatrix(self):
 1743         """
 1744         Subspace -> Field(Square)Matrix
 1745         """
 1746         if self.row == self.column:
 1747             self.__class__ = FieldSquareMatrix
 1748         else:
 1749             self.__class__ = FieldMatrix
 1750 
 1751     def isSubspace(self, other):
 1752         """ 
 1753         Check self is in other as subspace
 1754         """
 1755         try:
 1756             other.inverseImage(self)
 1757             return True
 1758         except:
 1759             return False
 1760 
 1761     def toBasis(self):
 1762         """
 1763         Change matrix to basis.
 1764         """
 1765         if not self.isbasis:
 1766             basis = self.image()
 1767             if not basis: # zero space
 1768                 basis = zeroMatrix(self.row, 1, self.coeff_ring)
 1769             self.compo = basis.compo
 1770             self.column = basis.column
 1771             self.isbasis = True
 1772 
 1773     def supplementBasis(self):     # Modified Algorithm 2.3.6 of Cohen's book
 1774         """
 1775         Return a basis of full space, which including self's column vectors.
 1776         """
 1777         self.toBasis()
 1778         if self.row < self.column:
 1779             raise MatrixSizeError()
 1780         n = self.row
 1781         k = self.column
 1782         M = self.copy()
 1783         pnt = range(1, self.row + 1)
 1784         for s in range(1, k + 1):
 1785             for t in range(s, n + 1):
 1786                 if M[t, s]:
 1787                     break
 1788             else: # zero space
 1789                 return unitMatrix(self.row, self.coeff_ring)
 1790             d = ring.inverse(M[t, s])
 1791             M[t, s] = 1
 1792             if t != s:
 1793                 pnt[t - 1] = pnt[s - 1]
 1794             for j in range(s + 1, k + 1):
 1795                 if t != s:
 1796                     tmp = M[s, j]
 1797                     M[s, j] = M[t, j]
 1798                     M[t, j] = tmp
 1799                 M[s, j] *= d
 1800                 for i in range(1, n + 1):
 1801                     if i != s and i != t:
 1802                         M[i, j] = M[i, j] - M[i, s] * M[s, j]
 1803         B = self.copy()
 1804         one = self.coeff_ring.one
 1805         zeros = [self.coeff_ring.zero] * n
 1806         for i in pnt[k: ]:
 1807             e_i = zeros
 1808             e_i[i - 1] = one
 1809             B.extendColumn(e_i)
 1810         return Subspace.fromMatrix(B, True)
 1811 
 1812     def sumOfSubspaces(self, other): # Algorithm 2.3.8 of Cohen's book
 1813         """
 1814         Return space which is sum of self and other.
 1815         """
 1816         if self.row != other.row:
 1817             raise MatrixSizeError()
 1818         N = self.copy()
 1819         N.extendColumn(other)
 1820         return Subspace.fromMatrix(N.image(), True)
 1821 
 1822     def intersectionOfSubspaces(self, other): # Algorithm 2.3.9 of Cohen's book
 1823         """
 1824         Return space which is intersection of self and other.
 1825         """
 1826         if self.row != other.row:
 1827             raise MatrixSizeError()
 1828         M1 = self.copy()
 1829         M1.extendColumn(other)
 1830         N = M1.kernel()
 1831         if not N:
 1832             zeroMatrix(self.row, 1, self.coeff_ring)
 1833         N1 = N.getBlock(1, 1, self.column, N.column) # N.column is dim(ker(M1))
 1834         return Subspace.fromMatrix((self * N1).image(), True)
 1835 
 1836 
 1837 # --------------------------------------------------------------------
 1838 #       the belows are not class methods
 1839 # --------------------------------------------------------------------
 1840 
 1841 def createMatrix(row, column=0, compo=0, coeff_ring=0):
 1842     """
 1843     generate new Matrix or SquareMatrix class.
 1844     """
 1845     if isinstance(compo, (ring.Ring, int)):
 1846         coeff_ring = compo
 1847         compo = 0
 1848     if isinstance(column, list):
 1849         compo = column
 1850         column = row
 1851     elif isinstance(column, ring.Ring):
 1852         coeff_ring = column
 1853         column = row
 1854     if coeff_ring != 0 and isinstance(coeff_ring, int):
 1855         coeff_ring = _getRing(coeff_ring)
 1856     if compo == 0:
 1857         return zeroMatrix(row, column, coeff_ring)
 1858     if coeff_ring == 0:
 1859         if isinstance(compo[0], (list, tuple)):
 1860             coeff_ring = ring.getRing(compo[0][0])
 1861         elif isinstance(compo[0], vector.Vector):
 1862             coeff_ring = ring.getRing(compo[0][1])
 1863         else:
 1864             coeff_ring = ring.getRing(compo[0])
 1865     if coeff_ring.isfield():
 1866         if row == column:
 1867             return FieldSquareMatrix(row, compo, coeff_ring)
 1868         else:
 1869             return FieldMatrix(row, column, compo, coeff_ring)
 1870     else:
 1871         if row == column:
 1872             return RingSquareMatrix(row, compo, coeff_ring)
 1873         else:
 1874             return RingMatrix(row, column, compo, coeff_ring)
 1875 
 1876 def unitMatrix(size, coeff=1):
 1877     """
 1878     return unit matrix of size.
 1879     coeff is subclass for ring.Ring or ring.Ring.one.
 1880     """
 1881     if isinstance(coeff, ring.Ring):
 1882         one = coeff.one
 1883         zero = coeff.zero
 1884     else:
 1885         one = coeff
 1886         coeff = ring.getRing(one)
 1887         zero = coeff.zero
 1888     unit_matrix = [one]
 1889     units = [zero] * size + [one]
 1890     for i in range(size - 1):
 1891         unit_matrix = unit_matrix + units
 1892     return createMatrix(size, size, unit_matrix, coeff)
 1893 
 1894 identityMatrix = unitMatrix
 1895 
 1896 def zeroMatrix(row, column=None, coeff=0):
 1897     """
 1898     return zero matrix.
 1899     coeff is subclass for ring.Ring or ring.Ring.zero.
 1900     """
 1901     if column == 0:
 1902         coeff = 0
 1903         column = row
 1904     if not(isinstance(column, (int, long))):
 1905         if column == None:
 1906             column = row
 1907         else:
 1908             coeff = column
 1909             column = row
 1910     if isinstance(coeff, ring.Ring):
 1911         zero = coeff.zero
 1912     else:
 1913         zero = coeff
 1914         coeff = ring.getRing(coeff)
 1915     zero_matrix = [zero] * (row * column)
 1916     return createMatrix(row, column, zero_matrix, coeff)
 1917 
 1918 
 1919 #--------------------------------------------------------------------
 1920 #   define exceptions
 1921 #--------------------------------------------------------------------
 1922 
 1923 class MatrixSizeError(Exception):
 1924     """Invalid input error for matrix size."""
 1925     pass
 1926 
 1927 class VectorsNotIndependent(Exception):
 1928     """Invalid input error because column vectors are linear dependent."""
 1929     pass
 1930 
 1931 class NoInverseImage(Exception):
 1932     """Invalid input error because self do not have inverse image for input."""
 1933     pass
 1934 
 1935 class NoInverse(Exception):
 1936     """Invalid input error because matrix is not invertible."""
 1937     pass
 1938 
 1939 ########## for utility
 1940 def _getRing(coeff_ring):
 1941     """
 1942     If 'coeff_ring' represents characteristic, return F_p or Z_n instance.
 1943     """
 1944     if isinstance(coeff_ring, int):
 1945         try:
 1946             import nzmath.finitefield
 1947             coeff_ring = nzmath.finitefield.FinitePrimeField(coeff_ring)
 1948         except:
 1949             import nzmath.intresidue
 1950             coeff_ring = nzmath.intresidue.IntegerResidueClassRing(coeff_ring)
 1951     return coeff_ring