X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fmatrix.cpp;h=f76aa9d768f0de53419c0285ebf9fc9485226596;hp=37716290f77317e4916b1705073ba3834332b1ce;hb=61e847cdfd8455a7ba67a0038add8dd35b4e0bff;hpb=5b090bb7e4951b48a28c20ba21bf9810c86eb0ca

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 37716290..f76aa9d7 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -26,9 +26,13 @@
 
 #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 {
@@ -77,9 +81,9 @@ const matrix & matrix::operator=(const matrix & other)
 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)
@@ -104,7 +108,7 @@ matrix::matrix(unsigned r, unsigned c)
     m.resize(r*c, _ex0());
 }
 
-// protected
+ // protected
 
 /** Ctor from representation, for internal use only. */
 matrix::matrix(unsigned r, unsigned c, const exvector & m2)
@@ -377,11 +381,13 @@ matrix matrix::mul(const matrix & other) const
         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);
@@ -466,9 +472,14 @@ ex matrix::determinant(void) const
         }
     }
     
-    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.
@@ -498,14 +509,9 @@ ex matrix::determinant(void) const
             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();
 }
 
 
@@ -533,24 +539,54 @@ ex matrix::trace(void) const
 }
 
 
-/** 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);
 }
 
 
@@ -633,7 +669,7 @@ ex matrix::ffe_get(unsigned r, unsigned c) const
 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"));
     
@@ -840,6 +876,7 @@ ex matrix::determinant_numeric(void) const
     ex det = _ex1();
     ex piv;
     
+    // standard Gauss method:
     for (unsigned r1=0; r1<row; ++r1) {
         int indx = tmp.pivot(r1);
         if (indx == -1)
@@ -859,67 +896,6 @@ ex matrix::determinant_numeric(void) const
 }
 
 
-/*  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")
@@ -930,7 +906,7 @@ ex matrix::determinant_minor_sparse(void) const
  *
  *  @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)
@@ -1039,45 +1015,110 @@ ex matrix::determinant_minor_dense(void) const
     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];
 }
 
 
@@ -1142,10 +1183,19 @@ int matrix::division_free_elimination(void)
  *  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)
@@ -1153,10 +1203,22 @@ int matrix::fraction_free_elimination(void)
         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();
         }
@@ -1213,6 +1275,29 @@ int matrix::pivot(unsigned ro, bool symbolic)
     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
 //////////