* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
+#include <iostream>
#include <algorithm>
#include <map>
#include <stdexcept>
#include "print.h"
#include "archive.h"
#include "utils.h"
-#include "debugmsg.h"
namespace GiNaC {
/** Default ctor. Initializes to 1 x 1-dimensional zero-matrix. */
matrix::matrix() : inherited(TINFO_matrix), row(1), col(1)
{
- debugmsg("matrix default ctor",LOGLEVEL_CONSTRUCT);
- m.push_back(_ex0());
+ m.push_back(_ex0);
}
void matrix::copy(const matrix & other)
matrix::matrix(unsigned r, unsigned c)
: inherited(TINFO_matrix), row(r), col(c)
{
- debugmsg("matrix ctor from unsigned,unsigned",LOGLEVEL_CONSTRUCT);
- m.resize(r*c, _ex0());
+ m.resize(r*c, _ex0);
}
// protected
/** Ctor from representation, for internal use only. */
matrix::matrix(unsigned r, unsigned c, const exvector & m2)
- : inherited(TINFO_matrix), row(r), col(c), m(m2)
-{
- debugmsg("matrix ctor from unsigned,unsigned,exvector",LOGLEVEL_CONSTRUCT);
-}
+ : inherited(TINFO_matrix), row(r), col(c), m(m2) {}
/** Construct matrix from (flat) list of elements. If the list has fewer
* elements than the matrix, the remaining matrix elements are set to zero.
matrix::matrix(unsigned r, unsigned c, const lst & l)
: inherited(TINFO_matrix), row(r), col(c)
{
- debugmsg("matrix ctor from unsigned,unsigned,lst",LOGLEVEL_CONSTRUCT);
- m.resize(r*c, _ex0());
+ m.resize(r*c, _ex0);
for (unsigned i=0; i<l.nops(); i++) {
unsigned x = i % c;
matrix::matrix(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
{
- debugmsg("matrix ctor from archive_node", LOGLEVEL_CONSTRUCT);
if (!(n.find_unsigned("row", row)) || !(n.find_unsigned("col", col)))
throw (std::runtime_error("unknown matrix dimensions in archive"));
m.reserve(row * col);
DEFAULT_UNARCHIVE(matrix)
//////////
-// functions overriding virtual functions from bases classes
+// functions overriding virtual functions from base classes
//////////
// public
void matrix::print(const print_context & c, unsigned level) const
{
- debugmsg("matrix print", LOGLEVEL_PRINT);
-
- if (is_of_type(c, print_tree)) {
+ if (is_a<print_tree>(c)) {
inherited::print(c, level);
} else {
+ if (is_a<print_python_repr>(c))
+ c.s << class_name() << '(';
+
c.s << "[";
for (unsigned y=0; y<row-1; ++y) {
c.s << "[";
m[row*col-1].print(c);
c.s << "]]";
+ if (is_a<print_python_repr>(c))
+ c.s << ')';
+
}
}
/** Evaluate matrix entry by entry. */
ex matrix::eval(int level) const
{
- debugmsg("matrix eval",LOGLEVEL_MEMBER_FUNCTION);
-
// check if we have to do anything at all
if ((level==1)&&(flags & status_flags::evaluated))
return *this;
for (unsigned c=0; c<col; ++c)
m2[r*col+c] = m[r*col+c].subs(ls, lr, no_pattern);
- return ex(matrix(row, col, m2)).bp->basic::subs(ls, lr, no_pattern);
+ return matrix(row, col, m2).basic::subs(ls, lr, no_pattern);
}
// protected
int matrix::compare_same_type(const basic & other) const
{
- GINAC_ASSERT(is_exactly_of_type(other, matrix));
- const matrix & o = static_cast<const matrix &>(other);
+ GINAC_ASSERT(is_exactly_a<matrix>(other));
+ const matrix &o = static_cast<const matrix &>(other);
// compare number of rows
if (row != o.rows())
bool matrix::match_same_type(const basic & other) const
{
- GINAC_ASSERT(is_exactly_of_type(other, matrix));
+ GINAC_ASSERT(is_exactly_a<matrix>(other));
const matrix & o = static_cast<const matrix &>(other);
// The number of rows and columns must be the same. This is necessary to
/** Automatic symbolic evaluation of an indexed matrix. */
ex matrix::eval_indexed(const basic & i) const
{
- GINAC_ASSERT(is_of_type(i, indexed));
- GINAC_ASSERT(is_ex_of_type(i.op(0), matrix));
+ GINAC_ASSERT(is_a<indexed>(i));
+ GINAC_ASSERT(is_a<matrix>(i.op(0)));
bool all_indices_unsigned = static_cast<const indexed &>(i).all_index_values_are(info_flags::nonnegint);
/** Sum of two indexed matrices. */
ex matrix::add_indexed(const ex & self, const ex & other) const
{
- GINAC_ASSERT(is_ex_of_type(self, indexed));
- GINAC_ASSERT(is_ex_of_type(self.op(0), matrix));
- GINAC_ASSERT(is_ex_of_type(other, indexed));
+ GINAC_ASSERT(is_a<indexed>(self));
+ GINAC_ASSERT(is_a<matrix>(self.op(0)));
+ GINAC_ASSERT(is_a<indexed>(other));
GINAC_ASSERT(self.nops() == 2 || self.nops() == 3);
// Only add two matrices
/** Product of an indexed matrix with a number. */
ex matrix::scalar_mul_indexed(const ex & self, const numeric & other) const
{
- GINAC_ASSERT(is_ex_of_type(self, indexed));
- GINAC_ASSERT(is_ex_of_type(self.op(0), matrix));
+ GINAC_ASSERT(is_a<indexed>(self));
+ GINAC_ASSERT(is_a<matrix>(self.op(0)));
GINAC_ASSERT(self.nops() == 2 || self.nops() == 3);
const matrix &self_matrix = ex_to<matrix>(self.op(0));
/** Contraction of an indexed matrix with something else. */
bool matrix::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
{
- GINAC_ASSERT(is_ex_of_type(*self, indexed));
- GINAC_ASSERT(is_ex_of_type(*other, indexed));
+ GINAC_ASSERT(is_a<indexed>(*self));
+ GINAC_ASSERT(is_a<indexed>(*other));
GINAC_ASSERT(self->nops() == 2 || self->nops() == 3);
- GINAC_ASSERT(is_ex_of_type(self->op(0), matrix));
+ GINAC_ASSERT(is_a<matrix>(self->op(0)));
// Only contract with other matrices
if (!is_ex_of_type(other->op(0), matrix))
*self = self_matrix.mul(other_matrix.transpose())(0, 0);
}
}
- *other = _ex1();
+ *other = _ex1;
return true;
} else { // vector * matrix
*self = indexed(self_matrix.mul(other_matrix), other->op(2));
else
*self = indexed(self_matrix.transpose().mul(other_matrix), other->op(2));
- *other = _ex1();
+ *other = _ex1;
return true;
}
*self = indexed(other_matrix.mul(self_matrix), other->op(1));
else
*self = indexed(other_matrix.mul(self_matrix.transpose()), other->op(1));
- *other = _ex1();
+ *other = _ex1;
return true;
}
}
// A_ij * B_jk = (A*B)_ik
if (is_dummy_pair(self->op(2), other->op(1))) {
*self = indexed(self_matrix.mul(other_matrix), self->op(1), other->op(2));
- *other = _ex1();
+ *other = _ex1;
return true;
}
// A_ij * B_kj = (A*Btrans)_ik
if (is_dummy_pair(self->op(2), other->op(2))) {
*self = indexed(self_matrix.mul(other_matrix.transpose()), self->op(1), other->op(1));
- *other = _ex1();
+ *other = _ex1;
return true;
}
// A_ji * B_jk = (Atrans*B)_ik
if (is_dummy_pair(self->op(1), other->op(1))) {
*self = indexed(self_matrix.transpose().mul(other_matrix), self->op(2), other->op(2));
- *other = _ex1();
+ *other = _ex1;
return true;
}
// A_ji * B_kj = (B*A)_ki
if (is_dummy_pair(self->op(1), other->op(2))) {
*self = indexed(other_matrix.mul(self_matrix), other->op(1), self->op(2));
- *other = _ex1();
+ *other = _ex1;
return true;
}
}
// 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"));
int sign;
sign = tmp.division_free_elimination(true);
if (sign==0)
- return _ex0();
+ return _ex0;
ex det = tmp.m[row*col-1];
// factor out accumulated bogus slag
for (unsigned d=0; d<row-2; ++d)
default: {
// This is the minor expansion scheme. We always develop such
// that the smallest minors (i.e, the trivial 1x1 ones) are on the
- // rightmost column. For this to be efficient it turns out that
- // the emptiest columns (i.e. the ones with most zeros) should be
- // the ones on the right hand side. Therefore we presort the
- // columns of the matrix:
+ // rightmost column. For this to be efficient, empirical tests
+ // have shown that the emptiest columns (i.e. the ones with most
+ // zeros) should be the ones on the right hand side -- although
+ // this might seem counter-intuitive (and in contradiction to some
+ // literature like the FORM manual). Please go ahead and test it
+ // if you don't believe me! Therefore we presort the columns of
+ // the matrix:
typedef std::pair<unsigned,unsigned> uintpair;
std::vector<uintpair> c_zeros; // number of zeros in column
for (unsigned c=0; c<col; ++c) {
++acc;
c_zeros.push_back(uintpair(acc,c));
}
- sort(c_zeros.begin(),c_zeros.end());
+ std::sort(c_zeros.begin(),c_zeros.end());
std::vector<unsigned> pre_sort;
for (std::vector<uintpair>::const_iterator i=c_zeros.begin(); i!=c_zeros.end(); ++i)
pre_sort.push_back(i->second);
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;
// First populate the identity matrix supposed to become the right hand side.
matrix identity(row,col);
for (unsigned i=0; i<row; ++i)
- identity(i,i) = _ex1();
+ identity(i,i) = _ex1;
// Populate a dummy matrix of variables, just because of compatibility with
// matrix::solve() which wants this (for compatibility with under-determined
// 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;
Pkey.push_back(i);
unsigned fc = 0; // controls logic for our strange flipper counter
do {
- det = _ex0();
+ det = _ex0;
for (unsigned r=0; r<n-c; ++r) {
// maybe there is nothing to do?
if (m[Pkey[r]*n+c].is_zero())
}
// fill up left hand side with zeros
for (unsigned c=0; c<=r1; ++c)
- this->m[r2*n+c] = _ex0();
+ this->m[r2*n+c] = _ex0;
}
if (det) {
// save space by deleting no longer needed elements
for (unsigned c=r0+1; c<n; ++c)
- this->m[r0*n+c] = _ex0();
+ this->m[r0*n+c] = _ex0;
}
++r0;
}
this->m[r2*n+c] = (this->m[r0*n+r1]*this->m[r2*n+c] - this->m[r2*n+r1]*this->m[r0*n+c]).expand();
// fill up left hand side with zeros
for (unsigned c=0; c<=r1; ++c)
- this->m[r2*n+c] = _ex0();
+ this->m[r2*n+c] = _ex0;
}
if (det) {
// save space by deleting no longer needed elements
for (unsigned c=r0+1; c<n; ++c)
- this->m[r0*n+c] = _ex0();
+ this->m[r0*n+c] = _ex0;
}
++r0;
}
}
// fill up left hand side with zeros
for (unsigned c=0; c<=r1; ++c)
- tmp_n.m[r2*n+c] = _ex0();
+ tmp_n.m[r2*n+c] = _ex0;
}
if ((r1<n-1)&&(r0<m-1)) {
// compute next iteration's divisor
if (det) {
// save space by deleting no longer needed elements
for (unsigned c=0; c<n; ++c) {
- tmp_n.m[r0*n+c] = _ex0();
- tmp_d.m[r0*n+c] = _ex1();
+ tmp_n.m[r0*n+c] = _ex0;
+ tmp_d.m[r0*n+c] = _ex1;
}
}
}
++k;
} else {
// search largest element in column co beginning at row ro
- GINAC_ASSERT(is_ex_of_type(this->m[k*col+co],numeric));
+ GINAC_ASSERT(is_a<numeric>(this->m[k*col+co]));
unsigned kmax = k+1;
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));
+ GINAC_ASSERT(is_a<numeric>(this->m[kmax*col+co]));
numeric tmp = ex_to<numeric>(this->m[kmax*col+co]);
if (abs(tmp) > mmax) {
mmax = tmp;
if (l.op(i).nops() > j)
m(i, j) = l.op(i).op(j);
else
- m(i, j) = _ex0();
+ m(i, j) = _ex0;
return m;
}
git://www.ginac.de / ginac.git / blobdiff