"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
565         """
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         """
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
1035         """
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         """
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
```