- rewrote binary scanning in matrix::pow(), stealing from CLN's expt_pos().

[ginac.git] / ginac / matrix.cpp

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 91382c2fd303612a1899a8f59c7d59b25b478880..65dfd4dc601e9c228f7329d6e5e18f6c471c0718 100644 (file)
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -141,7 +141,7 @@ void matrix::archive(archive_node &n) const
 DEFAULT_UNARCHIVE(matrix)
 
 //////////
-// functions overriding virtual functions from bases classes
+// functions overriding virtual functions from base classes
 //////////
 
 // public
@@ -415,10 +415,8 @@ bool matrix::contract_with(exvector::iterator self, exvector::iterator other, ex
        const matrix &other_matrix = ex_to<matrix>(other->op(0));
 
        if (self->nops() == 2) {
-               unsigned self_dim = (self_matrix.col == 1) ? self_matrix.row : self_matrix.col;
 
                if (other->nops() == 2) { // vector * vector (scalar product)
-                       unsigned other_dim = (other_matrix.col == 1) ? other_matrix.row : other_matrix.col;
 
                        if (self_matrix.col == 1) {
                                if (other_matrix.col == 1) {
@@ -601,32 +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 and
-                       // multiplies the factors whenever needed:
-                       while (b.compare(k)<=0) {
-                               b *= numeric(2);
-                               numeric r(mod(k,b));
-                               if (!r.is_zero()) {
-                                       k -= r;
-                                       result = result.mul(prod);
+                               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!=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"));