X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Fmatrix.cpp;h=60869f09fa709bbd17b051206fa4b32c29db0e73;hb=55ef6d50c8a068080d8fffb22e04358cdfd289b8;hp=ca574cc32fd9f7355fe378dded9c9bde50cf9da8;hpb=1b5601580eba6925b0af9408c56d6b2027b56013;p=ginac.git

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index ca574cc3..60869f09 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -599,35 +599,31 @@ matrix matrix::pow(const ex & expn) const
 		// Integer cases are computed by successive multiplication, using the
 		// obvious shortcut of storing temporaries, like A^4 == (A*A)*(A*A).
 		if (expn.info(info_flags::integer)) {
-			numeric k;
-			matrix prod(row,col);
+			numeric b = ex_to<numeric>(expn);
+			matrix A(row,col);
 			if (expn.info(info_flags::negative)) {
-				k = -ex_to<numeric>(expn);
-				prod = this->inverse();
+				b *= -1;
+				A = this->inverse();
 			} else {
-				k = ex_to<numeric>(expn);
-				prod = *this;
+				A = *this;
 			}
-			matrix result(row,col);
+			matrix C(row,col);
 			for (unsigned r=0; r<row; ++r)
-				result(r,r) = _ex1();
-			numeric b(1);
-			// This loop computes the representation of k in base 2 from right
-			// to left(!) and multiplies the factors whenever needed.  Note
+				C(r,r) = _ex1();
+			// 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.compare(k)<=0) {
-				b *= numeric(2);
-				numeric r(mod(k,b));
-				if (!r.is_zero()) {
-					k -= r;
-					result = result.mul(prod);
+			while (b!=1) {
+				if (b.is_odd()) {
+					C = C.mul(A);
+					b -= 1;
 				}
-				if (b.compare(k)<=0)
-					prod = prod.mul(prod);
+				b *= _num1_2();  // b /= 2, still integer.
+				A = A.mul(A);
 			}
-			return result;
+			return A.mul(C);
 		}
 	}
 	throw (std::runtime_error("matrix::pow(): don't know how to handle exponent"));
@@ -854,7 +850,7 @@ ex matrix::charpoly(const symbol & lambda) const
 	
 	bool numeric_flag = true;
 	exvector::const_iterator r = m.begin(), rend = m.end();
-	while (r != rend) {
+	while (r!=rend && numeric_flag==true) {
 		if (!r->info(info_flags::numeric))
 			numeric_flag = false;
 		++r;
@@ -965,7 +961,7 @@ matrix matrix::solve(const matrix & vars,
 	// Gather some statistical information about the augmented matrix:
 	bool numeric_flag = true;
 	exvector::const_iterator r = aug.m.begin(), rend = aug.m.end();
-	while (r != rend) {
+	while (r!=rend && numeric_flag==true) {
 		if (!r->info(info_flags::numeric))
 			numeric_flag = false;
 		++r;