#include "matrix.h"
#include "archive.h"
+#include "numeric.h"
+#include "lst.h"
#include "utils.h"
#include "debugmsg.h"
-#include "numeric.h"
+#include "power.h"
+#include "symbol.h"
+#include "normal.h"
#ifndef NO_NAMESPACE_GINAC
namespace GiNaC {
void matrix::copy(const matrix & other)
{
inherited::copy(other);
- row=other.row;
- col=other.col;
- m=other.m; // use STL's vector copying
+ row = other.row;
+ col = other.col;
+ m = other.m; // STL's vector copying invoked here
}
void matrix::destroy(bool call_parent)
m.resize(r*c, _ex0());
}
-// protected
+ // protected
/** Ctor from representation, for internal use only. */
matrix::matrix(unsigned r, unsigned c, const exvector & m2)
throw (std::logic_error("matrix::mul(): incompatible matrices"));
exvector prod(row*other.col);
- for (unsigned i=0; i<row; ++i) {
- for (unsigned j=0; j<other.col; ++j) {
- for (unsigned l=0; l<col; ++l) {
- prod[i*other.col+j] += m[i*col+l] * other.m[l*other.col+j];
- }
+
+ for (unsigned r1=0; r1<row; ++r1) {
+ for (unsigned c=0; c<col; ++c) {
+ if (m[r1*col+c].is_zero())
+ continue;
+ for (unsigned r2=0; r2<other.col; ++r2)
+ prod[r1*other.col+r2] += m[r1*col+c] * other.m[c*other.col+r2];
}
}
return matrix(row, other.col, prod);
}
}
- if (numeric_flag)
+ if (numeric_flag) // purely numeric matrix
return determinant_numeric();
+ // Does anybody really know when a matrix is sparse?
+ if (4*sparse_count<row*col) { // < row/2 nonzero elements average in a row
+ return determinant_bareiss(normal_flag);
+ }
+
// Now come the minor expansion schemes. 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.
result[r*col+c] = m[r*col+(*i)];
}
- if (0*sparse_count>row*col) { // MAGIC, maybe 10 some day?
- if (normal_flag)
- return sign*matrix(row,col,result).determinant_minor_sparse().normal();
- return sign*matrix(row,col,result).determinant_minor_sparse();
- }
if (normal_flag)
- return sign*matrix(row,col,result).determinant_minor_dense().normal();
- return sign*matrix(row,col,result).determinant_minor_dense();
+ return sign*matrix(row,col,result).determinant_minor().normal();
+ return sign*matrix(row,col,result).determinant_minor();
}
}
-/** Characteristic Polynomial. The characteristic polynomial of a matrix M is
- * defined as the determiant of (M - lambda * 1) where 1 stands for the unit
- * matrix of the same dimension as M. This method returns the characteristic
- * polynomial as a new expression.
+/** Characteristic Polynomial. Following mathematica notation the
+ * characteristic polynomial of a matrix M is defined as the determiant of
+ * (M - lambda * 1) where 1 stands for the unit matrix of the same dimension
+ * as M. Note that some CASs define it with a sign inside the determinant
+ * which gives rise to an overall sign if the dimension is odd. This method
+ * returns the characteristic polynomial collected in powers of lambda as a
+ * new expression.
*
* @return characteristic polynomial as new expression
* @exception logic_error (matrix not square)
* @see matrix::determinant() */
-ex matrix::charpoly(const ex & lambda) const
+ex matrix::charpoly(const symbol & lambda) const
{
if (row != col)
throw (std::logic_error("matrix::charpoly(): matrix not square"));
+ bool numeric_flag = true;
+ for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r) {
+ if (!(*r).info(info_flags::numeric)) {
+ numeric_flag = false;
+ }
+ }
+
+ // The pure numeric case is traditionally rather common. Hence, it is
+ // trapped and we use Leverrier's algorithm which goes as row^3 for
+ // every coefficient. The expensive section is the matrix multiplication.
+ if (numeric_flag) {
+ matrix B(*this);
+ ex c = B.trace();
+ ex poly = power(lambda,row)-c*power(lambda,row-1);
+ for (unsigned i=1; i<row; ++i) {
+ for (unsigned j=0; j<row; ++j)
+ B.m[j*col+j] -= c;
+ B = this->mul(B);
+ c = B.trace()/ex(i+1);
+ poly -= c*power(lambda,row-i-1);
+ }
+ if (row%2)
+ return -poly;
+ else
+ return poly;
+ }
+
matrix M(*this);
for (unsigned r=0; r<col; ++r)
M.m[r*col+r] -= lambda;
- return (M.determinant());
+ return M.determinant().collect(lambda);
}
matrix matrix::fraction_free_elim(const matrix & vars,
const matrix & rhs) const
{
- // FIXME: implement a Sasaki-Murao scheme which avoids division at all!
+ // FIXME: use implementation of matrix::fraction_free_elim
if ((row != rhs.row) || (col != vars.row) || (rhs.col != vars.col))
throw (std::logic_error("matrix::fraction_free_elim(): incompatible matrices"));
ex det = _ex1();
ex piv;
+ // standard Gauss method:
for (unsigned r1=0; r1<row; ++r1) {
int indx = tmp.pivot(r1);
if (indx == -1)
}
-/* Leverrier algorithm for large matrices having at least one symbolic entry.
- * This routine is only called internally by matrix::determinant(). The
- * algorithm is very bad for symbolic matrices since it returns expressions
- * that are quite hard to expand. */
-/*ex matrix::determinant_leverrier(const matrix & M)
- *{
- * GINAC_ASSERT(M.rows()==M.cols()); // cannot happen, just in case...
- *
- * matrix B(M);
- * matrix I(M.row, M.col);
- * ex c=B.trace();
- * for (unsigned i=1; i<M.row; ++i) {
- * for (unsigned j=0; j<M.row; ++j)
- * I.m[j*M.col+j] = c;
- * B = M.mul(B.sub(I));
- * c = B.trace()/ex(i+1);
- * }
- * if (M.row%2) {
- * return c;
- * } else {
- * return -c;
- * }
- *}*/
-
-
-ex matrix::determinant_minor_sparse(void) const
-{
- // for small matrices the algorithm does not make any sense:
- if (this->row==1)
- return m[0];
- if (this->row==2)
- return (m[0]*m[3]-m[2]*m[1]).expand();
- if (this->row==3)
- return (m[0]*m[4]*m[8]-m[0]*m[5]*m[7]-
- m[1]*m[3]*m[8]+m[2]*m[3]*m[7]+
- m[1]*m[5]*m[6]-m[2]*m[4]*m[6]).expand();
-
- ex det;
- matrix minorM(this->row-1,this->col-1);
- for (unsigned r1=0; r1<this->row; ++r1) {
- // shortcut if element(r1,0) vanishes
- if (m[r1*col].is_zero())
- continue;
- // assemble the minor matrix
- for (unsigned r=0; r<minorM.rows(); ++r) {
- for (unsigned c=0; c<minorM.cols(); ++c) {
- if (r<r1)
- minorM.set(r,c,m[r*col+c+1]);
- else
- minorM.set(r,c,m[(r+1)*col+c+1]);
- }
- }
- // recurse down and care for sign:
- if (r1%2)
- det -= m[r1*col] * minorM.determinant_minor_sparse();
- else
- det += m[r1*col] * minorM.determinant_minor_sparse();
- }
- return det.expand();
-}
-
/** Recursive determinant for small matrices having at least one symbolic
* entry. The basic algorithm, known as Laplace-expansion, is enhanced by
* some bookkeeping to avoid calculation of the same submatrices ("minors")
*
* @return the determinant as a new expression (in expanded form)
* @see matrix::determinant() */
-ex matrix::determinant_minor_dense(void) const
+ex matrix::determinant_minor(void) const
{
// for small matrices the algorithm does not make any sense:
if (this->row==1)
return det;
}
-
-/* Determinant using a simple Bareiss elimination scheme. Suited for
- * sparse matrices.
- *
- * @return the determinant as a new expression (in expanded form)
- * @see matrix::determinant() */
-ex matrix::determinant_bareiss(void) const
+/** Helper function to divide rational functions, as needed in any Bareiss
+ * elimination scheme over a quotient field.
+ *
+ * @see divide() */
+bool rat_divide(const ex & a, const ex & b, ex & q, bool check_args = true)
{
- matrix M(*this);
- int sign = M.fraction_free_elimination();
- if (sign)
- return sign*M(row-1,col-1);
- else
- return _ex0();
+ q = _ex0();
+ if (b.is_zero())
+ throw(std::overflow_error("rat_divide(): division by zero"));
+ if (a.is_zero())
+ return true;
+ if (is_ex_exactly_of_type(b, numeric)) {
+ q = a / b;
+ return true;
+ } else if (is_ex_exactly_of_type(a, numeric))
+ return false;
+ ex a_n = a.numer();
+ ex a_d = a.denom();
+ ex b_n = b.numer();
+ ex b_d = b.denom();
+ ex n; // new numerator
+ ex d; // new denominator
+ bool check = true;
+ check &= divide(a_n, b_n, n, check_args);
+ check &= divide(a_d, b_d, d, check_args);
+ q = n/d;
+ return check;
}
-/** Determinant built by application of the full permutation group. This
+/** Determinant computed by using fraction free elimination. This
* routine is only called internally by matrix::determinant().
- * NOTE: it is probably inefficient in all cases and may be eliminated. */
-ex matrix::determinant_perm(void) const
+ *
+ * @param normalize may be set to false only in integral domains. */
+ex matrix::determinant_bareiss(bool normalize) const
{
- if (rows()==1) // speed things up
+ if (rows()==1)
return m[0];
- ex det;
- ex term;
- vector<unsigned> sigma(col);
- for (unsigned i=0; i<col; ++i)
- sigma[i]=i;
+ int sign = 1;
+ ex divisor = 1;
+ ex dividend;
- do {
- term = (*this)(sigma[0],0);
- for (unsigned i=1; i<col; ++i)
- term *= (*this)(sigma[i],i);
- det += permutation_sign(sigma)*term;
- } while (next_permutation(sigma.begin(), sigma.end()));
+ // we populate a tmp matrix to subsequently operate on, it should
+ // be normalized even though this algorithm doesn't need GCDs since
+ // the elements of *this might be unnormalized, which complicates
+ // things:
+ matrix tmp(*this);
+ exvector::const_iterator i = m.begin();
+ exvector::iterator ti = tmp.m.begin();
+ for (; i!= m.end(); ++i, ++ti) {
+ if (normalize)
+ (*ti) = (*i).normal();
+ else
+ (*ti) = (*i);
+ }
- return det;
+ for (unsigned r1=0; r1<row-1; ++r1) {
+ int indx = tmp.pivot(r1);
+ if (indx==-1)
+ return _ex0();
+ if (indx>0)
+ sign = -sign;
+ if (r1>0) {
+ divisor = tmp.m[(r1-1)*col+(r1-1)].expand();
+ // delete the elements we don't need anymore:
+ for (unsigned c=0; c<col; ++c)
+ tmp.m[(r1-1)*col+c] = _ex0();
+ }
+ for (unsigned r2=r1+1; r2<row; ++r2) {
+ for (unsigned c=r1+1; c<col; ++c) {
+ lst srl; // symbol replacement list for .to_rational()
+ dividend = (tmp.m[r1*tmp.col+r1]*tmp.m[r2*tmp.col+c]
+ -tmp.m[r2*tmp.col+r1]*tmp.m[r1*tmp.col+c]).expand();
+ if (normalize) {
+#ifdef DO_GINAC_ASSERT
+ GINAC_ASSERT(rat_divide(dividend.to_rational(srl),
+ divisor.to_rational(srl),
+ tmp.m[r2*tmp.col+c],true));
+#else
+ rat_divide(dividend.to_rational(srl),
+ divisor.to_rational(srl),
+ tmp.m[r2*tmp.col+c],false);
+#endif
+ }
+ else {
+#ifdef DO_GINAC_ASSERT
+ GINAC_ASSERT(divide(dividend.to_rational(srl),
+ divisor.to_rational(srl),
+ tmp.m[r2*tmp.col+c],true));
+#else
+ divide(dividend.to_rational(srl),
+ divisor.to_rational(srl),
+ tmp.m[r2*tmp.col+c],false);
+#endif
+ }
+ tmp.m[r2*tmp.col+c] = tmp.m[r2*tmp.col+c].subs(srl);
+ }
+ for (unsigned c=0; c<=r1; ++c)
+ tmp.m[r2*tmp.col+c] = _ex0();
+ }
+ }
+
+ return sign*tmp.m[tmp.row*tmp.col-1];
}
* number of rows was swapped and 0 if the matrix is singular. */
int matrix::fraction_free_elimination(void)
{
+ ensure_if_modifiable();
+
+ // first normal all elements:
+ for (exvector::iterator i=m.begin(); i!=m.end(); ++i)
+ (*i) = (*i).normal();
+
+ // FIXME: this is unfinished, once matrix::determinant_bareiss is
+ // bulletproof, some code ought to be copy from there to here.
int sign = 1;
ex divisor = 1;
+ ex dividend;
+ lst srl; // symbol replacement list for .to_rational()
- ensure_if_modifiable();
for (unsigned r1=0; r1<row-1; ++r1) {
int indx = pivot(r1);
if (indx==-1)
if (indx>0)
sign = -sign;
if (r1>0)
- divisor = this->m[(r1-1)*col + (r1-1)];
+ divisor = this->m[(r1-1)*col+(r1-1)].expand();
for (unsigned r2=r1+1; r2<row; ++r2) {
- for (unsigned c=r1+1; c<col; ++c)
- this->m[r2*col+c] = ((this->m[r1*col+r1]*this->m[r2*col+c] - this->m[r2*col+r1]*this->m[r1*col+c])/divisor).normal();
+ for (unsigned c=r1+1; c<col; ++c) {
+ dividend = (this->m[r1*col+r1]*this->m[r2*col+c]
+ -this->m[r2*col+r1]*this->m[r1*col+c]).expand();
+#ifdef DO_GINAC_ASSERT
+ GINAC_ASSERT(divide(dividend.to_rational(srl),
+ divisor.to_rational(srl),
+ this->m[r2*col+c]));
+#else
+ divide(dividend.to_rational(srl),
+ divisor.to_rational(srl),
+ this->m[r2*col+c]);
+#endif // DO_GINAC_ASSERT
+ this->m[r2*col+c] = this->m[r2*col+c].subs(srl);
+ }
for (unsigned c=0; c<=r1; ++c)
this->m[r2*col+c] = _ex0();
}
return 0;
}
+/** Convert list of lists to matrix. */
+ex lst_to_matrix(const ex &l)
+{
+ if (!is_ex_of_type(l, lst))
+ throw(std::invalid_argument("argument to lst_to_matrix() must be a lst"));
+
+ // Find number of rows and columns
+ unsigned rows = l.nops(), cols = 0, i, j;
+ for (i=0; i<rows; i++)
+ if (l.op(i).nops() > cols)
+ cols = l.op(i).nops();
+
+ // Allocate and fill matrix
+ matrix &m = *new matrix(rows, cols);
+ for (i=0; i<rows; i++)
+ for (j=0; j<cols; j++)
+ if (l.op(i).nops() > j)
+ m.set(i, j, l.op(i).op(j));
+ else
+ m.set(i, j, ex(0));
+ return m;
+}
+
//////////
// global constants
//////////
git://www.ginac.de / ginac.git / blobdiff