- comment: elaborate on matrix::pow().

[ginac.git] / ginac / matrix.cpp

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 1f5e9af50f939273a17727f6d7b78bd83682c233..ca574cc32fd9f7355fe378dded9c9bde50cf9da8 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
@@ -612,8 +612,11 @@ matrix matrix::pow(const ex & expn) const
                        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:
+                       // This loop computes the representation of k 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));