- Remove the -ansi compiler switch in the example since it doesn't

[ginac.git] / ginac / matrix.cpp

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index ca574cc32fd9f7355fe378dded9c9bde50cf9da8..60869f09fa709bbd17b051206fa4b32c29db0e73 100644 (file)
--- 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;