From: Richard Kreckel <Richard.Kreckel@uni-mainz.de>
Date: Fri, 8 Mar 2002 20:32:59 +0000 (+0000)
Subject: * matrix:;pow(): Fix exponent==0 case.
X-Git-Tag: release_1-0-7~2
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=b3f0d2416a7a0195446e410ebba10ec2c51cbd37

* matrix:;pow(): Fix exponent==0 case.
* matrix::print(): Implement print_latex.
* Add unit_matrix(unsigend) and unit_matrix(unsigned,unsigned).
---

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index f92a7c56..8daa11db 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -150,23 +150,38 @@ void matrix::print(const print_context & c, unsigned level) const
 		if (is_a<print_python_repr>(c))
 			c.s << class_name() << '(';
 
-		c.s << "[";
-		for (unsigned y=0; y<row-1; ++y) {
+		if (is_a<print_latex>(c))
+			c.s << "\\left(\\begin{array}{" << std::string(col,'c') << "}";
+		else
 			c.s << "[";
-			for (unsigned x=0; x<col-1; ++x) {
-				m[y*col+x].print(c);
-				c.s << ",";
+
+		for (unsigned ro=0; ro<row; ++ro) {
+			if (!is_a<print_latex>(c))
+				c.s << "[";
+			for (unsigned co=0; co<col; ++co) {
+				m[ro*col+co].print(c);
+				if (co<col-1) {
+					if (is_a<print_latex>(c))
+						c.s << "&";
+					else
+						c.s << ",";
+				} else {
+					if (!is_a<print_latex>(c))
+						c.s << "]";
+				}
+			}
+			if (ro<row-1) {
+				if (is_a<print_latex>(c))
+					c.s << "\\\\";
+				else
+					c.s << ",";
 			}
-			m[col*(y+1)-1].print(c);
-			c.s << "],";
-		}
-		c.s << "[";
-		for (unsigned x=0; x<col-1; ++x) {
-			m[(row-1)*col+x].print(c);
-			c.s << ",";
 		}
-		m[row*col-1].print(c);
-		c.s << "]]";
+
+		if (is_a<print_latex>(c))
+			c.s << "\\end{array}\\right)";
+		else
+			c.s << "]";
 
 		if (is_a<print_python_repr>(c))
 			c.s << ')';
@@ -605,17 +620,19 @@ matrix matrix::pow(const ex & expn) const
 			matrix C(row,col);
 			for (unsigned r=0; r<row; ++r)
 				C(r,r) = _ex1;
+			if (b.is_zero())
+				return C;
 			// This loop computes the representation of b in base 2 from right
 			// to left and multiplies the factors whenever needed.  Note
 			// that this is not entirely optimal but close to optimal and
 			// "better" algorithms are much harder to implement.  (See Knuth,
 			// TAoCP2, section "Evaluation of Powers" for a good discussion.)
-			while (b!=1) {
+			while (b!=_num1) {
 				if (b.is_odd()) {
 					C = C.mul(A);
-					b -= 1;
+					--b;
 				}
-				b *= _num1_2;  // b /= 2, still integer.
+				b /= _num2;  // still integer.
 				A = A.mul(A);
 			}
 			return A.mul(C);
@@ -1460,4 +1477,12 @@ ex diag_matrix(const lst & l)
 	return m;
 }
 
+ex unit_matrix(unsigned r, unsigned c)
+{
+	matrix Id(r,c);
+	for (unsigned i=0; i<r && i<c; ++i)
+		Id(i,i) = _ex1;
+	return Id;
+}
+
 } // namespace GiNaC
diff --git a/ginac/matrix.h b/ginac/matrix.h
index 26783eab..518684ca 100644
--- a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ -144,6 +144,13 @@ extern ex lst_to_matrix(const lst & l);
 /** Convert list of diagonal elements to matrix. */
 extern ex diag_matrix(const lst & l);
 
+/** Create a r times c unit matrix. */
+extern ex unit_matrix(unsigned r, unsigned c);
+
+/** Create a x times x unit matrix. */
+inline ex unit_matrix(unsigned x)
+{ return unit_matrix(x, x); }
+
 } // namespace GiNaC
 
 #endif // ndef __GINAC_MATRIX_H__