"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