* Avoid getrusage(2) on systems that don't have it (by ASheplyakov Alexei

[ginac.git] / ginac / matrix.cpp

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index e24f644741d6be1745c610020d314b5c338cf7ac..710190376a1c28b59c787d902655aa24d97ccf73 100644 (file)
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -3,7 +3,7 @@
  *  Implementation of symbolic matrices */
 
 /*
- *  GiNaC Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -17,7 +17,7 @@
  *
  *  You should have received a copy of the GNU General Public License
  *  along with this program; if not, write to the Free Software
- *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
 #include <string>
@@ -647,12 +647,12 @@ matrix matrix::pow(const ex & expn) const
                        // 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!=_num1) {
+                       while (b!=*_num1_p) {
                                if (b.is_odd()) {
                                        C = C.mul(A);
                                        --b;
                                }
-                               b /= _num2;  // still integer.
+                               b /= *_num2_p;  // still integer.
                                A = A.mul(A);
                        }
                        return A.mul(C);
@@ -730,12 +730,12 @@ ex matrix::determinant(unsigned algo) const
        unsigned sparse_count = 0;  // counts non-zero elements
        exvector::const_iterator r = m.begin(), rend = m.end();
        while (r != rend) {
+               if (!r->info(info_flags::numeric))
+                       numeric_flag = false;
                exmap srl;  // symbol replacement list
                ex rtest = r->to_rational(srl);
                if (!rtest.is_zero())
                        ++sparse_count;
-               if (!rtest.info(info_flags::numeric))
-                       numeric_flag = false;
                if (!rtest.info(info_flags::crational_polynomial) &&
                         rtest.info(info_flags::rational_function))
                        normal_flag = true;
@@ -1205,7 +1205,7 @@ ex matrix::determinant_minor() const
                                        Pkey[j] = Pkey[j-1]+1;
                } while(fc);
                // next column, so change the role of A and B:
-               A = B;
+               A.swap(B);
                B.clear();
        }
        
@@ -1598,4 +1598,52 @@ ex symbolic_matrix(unsigned r, unsigned c, const std::string & base_name, const
        return M;
 }
 
+ex reduced_matrix(const matrix& m, unsigned r, unsigned c)
+{
+       if (r+1>m.rows() || c+1>m.cols() || m.cols()<2 || m.rows()<2)
+               throw std::runtime_error("minor_matrix(): index out of bounds");
+
+       const unsigned rows = m.rows()-1;
+       const unsigned cols = m.cols()-1;
+       matrix &M = *new matrix(rows, cols);
+       M.setflag(status_flags::dynallocated | status_flags::evaluated);
+
+       unsigned ro = 0;
+       unsigned ro2 = 0;
+       while (ro2<rows) {
+               if (ro==r)
+                       ++ro;
+               unsigned co = 0;
+               unsigned co2 = 0;
+               while (co2<cols) {
+                       if (co==c)
+                               ++co;
+                       M(ro2,co2) = m(ro, co);
+                       ++co;
+                       ++co2;
+               }
+               ++ro;
+               ++ro2;
+       }
+
+       return M;
+}
+
+ex sub_matrix(const matrix&m, unsigned r, unsigned nr, unsigned c, unsigned nc)
+{
+       if (r+nr>m.rows() || c+nc>m.cols())
+               throw std::runtime_error("sub_matrix(): index out of bounds");
+
+       matrix &M = *new matrix(nr, nc);
+       M.setflag(status_flags::dynallocated | status_flags::evaluated);
+
+       for (unsigned ro=0; ro<nr; ++ro) {
+               for (unsigned co=0; co<nc; ++co) {
+                       M(ro,co) = m(ro+r,co+c);
+               }
+       }
+
+       return M;
+}
+
 } // namespace GiNaC