- The default implementations of evalf(), diff(), normal() and expand() use

[ginac.git] / ginac / matrix.cpp

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 9998a3918ab7d50eb3c4685aad29b952db04c392..9f9f67a82c83b62e76431de79ea4bf4054a05047 100644 (file)
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -198,16 +198,6 @@ ex & matrix::let_op(int i)
        return m[i];
 }
 
-/** expands the elements of a matrix entry by entry. */
-ex matrix::expand(unsigned options) const
-{
-       exvector tmp(row*col);
-       for (unsigned i=0; i<row*col; ++i)
-               tmp[i] = m[i].expand(options);
-       
-       return matrix(row, col, tmp);
-}
-
 /** Evaluate matrix entry by entry. */
 ex matrix::eval(int level) const
 {
@@ -232,30 +222,6 @@ ex matrix::eval(int level) const
                                                                                           status_flags::evaluated );
 }
 
-/** Evaluate matrix numerically entry by entry. */
-ex matrix::evalf(int level) const
-{
-       debugmsg("matrix evalf",LOGLEVEL_MEMBER_FUNCTION);
-               
-       // check if we have to do anything at all
-       if (level==1)
-               return *this;
-       
-       // emergency break
-       if (level == -max_recursion_level) {
-               throw (std::runtime_error("matrix::evalf(): recursion limit exceeded"));
-       }
-       
-       // evalf() entry by entry
-       exvector m2(row*col);
-       --level;
-       for (unsigned r=0; r<row; ++r)
-               for (unsigned c=0; c<col; ++c)
-                       m2[r*col+c] = m[r*col+c].evalf(level);
-       
-       return matrix(row, col, m2);
-}
-
 ex matrix::subs(const lst & ls, const lst & lr, bool no_pattern) const
 {
        exvector m2(row * col);
@@ -308,7 +274,7 @@ ex matrix::eval_indexed(const basic & i) const
                if (row != 1 && col != 1)
                        throw (std::runtime_error("matrix::eval_indexed(): vector must have exactly 1 index"));
 
-               const idx & i1 = ex_to_idx(i.op(1));
+               const idx & i1 = ex_to<idx>(i.op(1));
 
                if (col == 1) {
 
@@ -318,7 +284,7 @@ ex matrix::eval_indexed(const basic & i) const
 
                        // Index numeric -> return vector element
                        if (all_indices_unsigned) {
-                               unsigned n1 = ex_to_numeric(i1.get_value()).to_int();
+                               unsigned n1 = ex_to<numeric>(i1.get_value()).to_int();
                                if (n1 >= row)
                                        throw (std::runtime_error("matrix::eval_indexed(): value of index exceeds number of vector elements"));
                                return (*this)(n1, 0);
@@ -332,7 +298,7 @@ ex matrix::eval_indexed(const basic & i) const
 
                        // Index numeric -> return vector element
                        if (all_indices_unsigned) {
-                               unsigned n1 = ex_to_numeric(i1.get_value()).to_int();
+                               unsigned n1 = ex_to<numeric>(i1.get_value()).to_int();
                                if (n1 >= col)
                                        throw (std::runtime_error("matrix::eval_indexed(): value of index exceeds number of vector elements"));
                                return (*this)(0, n1);
@@ -342,8 +308,8 @@ ex matrix::eval_indexed(const basic & i) const
        } else if (i.nops() == 3) {
 
                // Two indices
-               const idx & i1 = ex_to_idx(i.op(1));
-               const idx & i2 = ex_to_idx(i.op(2));
+               const idx & i1 = ex_to<idx>(i.op(1));
+               const idx & i2 = ex_to<idx>(i.op(2));
 
                if (!i1.get_dim().is_equal(row))
                        throw (std::runtime_error("matrix::eval_indexed(): dimension of first index must match number of rows"));
@@ -356,7 +322,7 @@ ex matrix::eval_indexed(const basic & i) const
 
                // Both indices numeric -> return matrix element
                if (all_indices_unsigned) {
-                       unsigned n1 = ex_to_numeric(i1.get_value()).to_int(), n2 = ex_to_numeric(i2.get_value()).to_int();
+                       unsigned n1 = ex_to<numeric>(i1.get_value()).to_int(), n2 = ex_to<numeric>(i2.get_value()).to_int();
                        if (n1 >= row)
                                throw (std::runtime_error("matrix::eval_indexed(): value of first index exceeds number of rows"));
                        if (n2 >= col)
@@ -382,8 +348,8 @@ ex matrix::add_indexed(const ex & self, const ex & other) const
        if (is_ex_of_type(other.op(0), matrix)) {
                GINAC_ASSERT(other.nops() == 2 || other.nops() == 3);
 
-               const matrix &self_matrix = ex_to_matrix(self.op(0));
-               const matrix &other_matrix = ex_to_matrix(other.op(0));
+               const matrix &self_matrix = ex_to<matrix>(self.op(0));
+               const matrix &other_matrix = ex_to<matrix>(other.op(0));
 
                if (self.nops() == 2 && other.nops() == 2) { // vector + vector
 
@@ -413,7 +379,7 @@ ex matrix::scalar_mul_indexed(const ex & self, const numeric & other) const
        GINAC_ASSERT(is_ex_of_type(self.op(0), matrix));
        GINAC_ASSERT(self.nops() == 2 || self.nops() == 3);
 
-       const matrix &self_matrix = ex_to_matrix(self.op(0));
+       const matrix &self_matrix = ex_to<matrix>(self.op(0));
 
        if (self.nops() == 2)
                return indexed(self_matrix.mul(other), self.op(1));
@@ -435,8 +401,8 @@ bool matrix::contract_with(exvector::iterator self, exvector::iterator other, ex
 
        GINAC_ASSERT(other->nops() == 2 || other->nops() == 3);
 
-       const matrix &self_matrix = ex_to_matrix(self->op(0));
-       const matrix &other_matrix = ex_to_matrix(other->op(0));
+       const matrix &self_matrix = ex_to<matrix>(self->op(0));
+       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;
@@ -628,18 +594,18 @@ matrix matrix::pow(const ex & expn) const
                        numeric k;
                        matrix prod(row,col);
                        if (expn.info(info_flags::negative)) {
-                               k = -ex_to_numeric(expn);
+                               k = -ex_to<numeric>(expn);
                                prod = this->inverse();
                        } else {
-                               k = ex_to_numeric(expn);
+                               k = ex_to<numeric>(expn);
                                prod = *this;
                        }
                        matrix result(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:
+                       // 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));
@@ -647,7 +613,8 @@ matrix matrix::pow(const ex & expn) const
                                        k -= r;
                                        result = result.mul(prod);
                                }
-                               prod = prod.mul(prod);
+                               if (b.compare(k)<=0)
+                                       prod = prod.mul(prod);
                        }
                        return result;
                }
@@ -681,7 +648,6 @@ ex & matrix::operator() (unsigned ro, unsigned co)
                throw (std::range_error("matrix::operator(): index out of range"));
 
        ensure_if_modifiable();
-       clearflag(status_flags::hash_calculated);
        return m[ro*col+co];
 }
 
@@ -1425,10 +1391,10 @@ int matrix::pivot(unsigned ro, unsigned co, bool symbolic)
                // search largest element in column co beginning at row ro
                GINAC_ASSERT(is_ex_of_type(this->m[k*col+co],numeric));
                unsigned kmax = k+1;
-               numeric mmax = abs(ex_to_numeric(m[kmax*col+co]));
+               numeric mmax = abs(ex_to<numeric>(m[kmax*col+co]));
                while (kmax<row) {
                        GINAC_ASSERT(is_ex_of_type(this->m[kmax*col+co],numeric));
-                       numeric tmp = ex_to_numeric(this->m[kmax*col+co]);
+                       numeric tmp = ex_to<numeric>(this->m[kmax*col+co]);
                        if (abs(tmp) > mmax) {
                                mmax = tmp;
                                k = kmax;