+/** Perform the steps of an ordinary Gaussian elimination to bring the m x n
+ * matrix into an upper echelon form. The algorithm is ok for matrices
+ * with numeric coefficients but quite unsuited for symbolic matrices.
+ *
+ * @param det may be set to true to save a lot of space if one is only
+ * interested in the diagonal elements (i.e. for calculating determinants).
+ * The others are set to zero in this case.
+ * @return sign is 1 if an even number of rows was swapped, -1 if an odd
+ * number of rows was swapped and 0 if the matrix is singular. */
+int matrix::gauss_elimination(const bool det)
+{
+ ensure_if_modifiable();
+ const unsigned m = this->rows();
+ const unsigned n = this->cols();
+ GINAC_ASSERT(!det || n==m);
+ int sign = 1;
+
+ unsigned r0 = 0;
+ for (unsigned r1=0; (r1<n-1)&&(r0<m-1); ++r1) {
+ int indx = pivot(r0, r1, true);
+ if (indx == -1) {
+ sign = 0;
+ if (det)
+ return 0; // leaves *this in a messy state
+ }
+ if (indx>=0) {
+ if (indx > 0)
+ sign = -sign;
+ for (unsigned r2=r0+1; r2<m; ++r2) {
+ if (!this->m[r2*n+r1].is_zero()) {
+ // yes, there is something to do in this row
+ ex piv = this->m[r2*n+r1] / this->m[r0*n+r1];
+ for (unsigned c=r1+1; c<n; ++c) {
+ this->m[r2*n+c] -= piv * this->m[r0*n+c];
+ if (!this->m[r2*n+c].info(info_flags::numeric))
+ this->m[r2*n+c] = this->m[r2*n+c].normal();
+ }
+ }
+ // fill up left hand side with zeros
+ for (unsigned c=0; c<=r1; ++c)
+ 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;
+ }
+ ++r0;
+ }
+ }
+
+ return sign;
+}
+
+
+/** Perform the steps of division free elimination to bring the m x n matrix
+ * into an upper echelon form.
+ *
+ * @param det may be set to true to save a lot of space if one is only
+ * interested in the diagonal elements (i.e. for calculating determinants).
+ * The others are set to zero in this case.
+ * @return sign is 1 if an even number of rows was swapped, -1 if an odd
+ * number of rows was swapped and 0 if the matrix is singular. */
+int matrix::division_free_elimination(const bool det)
+{
+ ensure_if_modifiable();
+ const unsigned m = this->rows();
+ const unsigned n = this->cols();
+ GINAC_ASSERT(!det || n==m);
+ int sign = 1;
+
+ unsigned r0 = 0;
+ for (unsigned r1=0; (r1<n-1)&&(r0<m-1); ++r1) {
+ int indx = pivot(r0, r1, true);
+ if (indx==-1) {
+ sign = 0;
+ if (det)
+ return 0; // leaves *this in a messy state
+ }
+ if (indx>=0) {
+ if (indx>0)
+ sign = -sign;
+ for (unsigned r2=r0+1; r2<m; ++r2) {
+ for (unsigned c=r1+1; c<n; ++c)
+ 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;
+ }
+ if (det) {
+ // save space by deleting no longer needed elements
+ for (unsigned c=r0+1; c<n; ++c)
+ this->m[r0*n+c] = _ex0;
+ }
+ ++r0;
+ }
+ }
+
+ return sign;
+}
+
+
+/** Perform the steps of Bareiss' one-step fraction free elimination to bring
+ * the matrix into an upper echelon form. Fraction free elimination means
+ * that divide is used straightforwardly, without computing GCDs first. This
+ * is possible, since we know the divisor at each step.
+ *
+ * @param det may be set to true to save a lot of space if one is only
+ * interested in the last element (i.e. for calculating determinants). The
+ * others are set to zero in this case.
+ * @return sign is 1 if an even number of rows was swapped, -1 if an odd
+ * number of rows was swapped and 0 if the matrix is singular. */
+int matrix::fraction_free_elimination(const bool det)
+{
+ // Method:
+ // (single-step fraction free elimination scheme, already known to Jordan)
+ //
+ // Usual division-free elimination sets m[0](r,c) = m(r,c) and then sets
+ // m[k+1](r,c) = m[k](k,k) * m[k](r,c) - m[k](r,k) * m[k](k,c).
+ //
+ // Bareiss (fraction-free) elimination in addition divides that element
+ // by m[k-1](k-1,k-1) for k>1, where it can be shown by means of the
+ // Sylvester identity that this really divides m[k+1](r,c).
+ //
+ // We also allow rational functions where the original prove still holds.
+ // However, we must care for numerator and denominator separately and
+ // "manually" work in the integral domains because of subtle cancellations
+ // (see below). This blows up the bookkeeping a bit and the formula has
+ // to be modified to expand like this (N{x} stands for numerator of x,
+ // D{x} for denominator of x):
+ // N{m[k+1](r,c)} = N{m[k](k,k)}*N{m[k](r,c)}*D{m[k](r,k)}*D{m[k](k,c)}
+ // -N{m[k](r,k)}*N{m[k](k,c)}*D{m[k](k,k)}*D{m[k](r,c)}
+ // D{m[k+1](r,c)} = D{m[k](k,k)}*D{m[k](r,c)}*D{m[k](r,k)}*D{m[k](k,c)}
+ // where for k>1 we now divide N{m[k+1](r,c)} by
+ // N{m[k-1](k-1,k-1)}
+ // and D{m[k+1](r,c)} by
+ // D{m[k-1](k-1,k-1)}.
+
+ ensure_if_modifiable();
+ const unsigned m = this->rows();
+ const unsigned n = this->cols();
+ GINAC_ASSERT(!det || n==m);
+ int sign = 1;
+ if (m==1)
+ return 1;
+ ex divisor_n = 1;
+ ex divisor_d = 1;
+ ex dividend_n;
+ ex dividend_d;
+
+ // We populate temporary matrices to subsequently operate on. There is
+ // one holding numerators and another holding denominators of entries.
+ // This is a must since the evaluator (or even earlier mul's constructor)
+ // might cancel some trivial element which causes divide() to fail. The
+ // elements are normalized first (yes, even though this algorithm doesn't
+ // need GCDs) since the elements of *this might be unnormalized, which
+ // makes things more complicated than they need to be.
+ matrix tmp_n(*this);
+ matrix tmp_d(m,n); // for denominators, if needed
+ lst srl; // symbol replacement list
+ exvector::const_iterator cit = this->m.begin(), citend = this->m.end();
+ exvector::iterator tmp_n_it = tmp_n.m.begin(), tmp_d_it = tmp_d.m.begin();
+ while (cit != citend) {
+ ex nd = cit->normal().to_rational(srl).numer_denom();
+ ++cit;
+ *tmp_n_it++ = nd.op(0);
+ *tmp_d_it++ = nd.op(1);
+ }
+
+ unsigned r0 = 0;
+ for (unsigned r1=0; (r1<n-1)&&(r0<m-1); ++r1) {
+ int indx = tmp_n.pivot(r0, r1, true);
+ if (indx==-1) {
+ sign = 0;
+ if (det)
+ return 0;
+ }
+ if (indx>=0) {
+ if (indx>0) {
+ sign = -sign;
+ // tmp_n's rows r0 and indx were swapped, do the same in tmp_d:
+ for (unsigned c=r1; c<n; ++c)
+ tmp_d.m[n*indx+c].swap(tmp_d.m[n*r0+c]);
+ }
+ for (unsigned r2=r0+1; r2<m; ++r2) {
+ for (unsigned c=r1+1; c<n; ++c) {
+ dividend_n = (tmp_n.m[r0*n+r1]*tmp_n.m[r2*n+c]*
+ tmp_d.m[r2*n+r1]*tmp_d.m[r0*n+c]
+ -tmp_n.m[r2*n+r1]*tmp_n.m[r0*n+c]*
+ tmp_d.m[r0*n+r1]*tmp_d.m[r2*n+c]).expand();
+ dividend_d = (tmp_d.m[r2*n+r1]*tmp_d.m[r0*n+c]*
+ tmp_d.m[r0*n+r1]*tmp_d.m[r2*n+c]).expand();
+ bool check = divide(dividend_n, divisor_n,
+ tmp_n.m[r2*n+c], true);
+ check &= divide(dividend_d, divisor_d,
+ tmp_d.m[r2*n+c], true);
+ GINAC_ASSERT(check);
+ }
+ // fill up left hand side with zeros
+ for (unsigned c=0; c<=r1; ++c)
+ tmp_n.m[r2*n+c] = _ex0;
+ }
+ if ((r1<n-1)&&(r0<m-1)) {
+ // compute next iteration's divisor
+ divisor_n = tmp_n.m[r0*n+r1].expand();
+ divisor_d = tmp_d.m[r0*n+r1].expand();
+ 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;
+ }
+ }
+ }
+ ++r0;
+ }
+ }
+ // repopulate *this matrix:
+ exvector::iterator it = this->m.begin(), itend = this->m.end();
+ tmp_n_it = tmp_n.m.begin();
+ tmp_d_it = tmp_d.m.begin();
+ while (it != itend)
+ *it++ = ((*tmp_n_it++)/(*tmp_d_it++)).subs(srl, subs_options::no_pattern);
+
+ return sign;
+}
+
+
+/** Partial pivoting method for matrix elimination schemes.
+ * Usual pivoting (symbolic==false) returns the index to the element with the
+ * largest absolute value in column ro and swaps the current row with the one
+ * where the element was found. With (symbolic==true) it does the same thing
+ * with the first non-zero element.
+ *
+ * @param ro is the row from where to begin
+ * @param co is the column to be inspected
+ * @param symbolic signal if we want the first non-zero element to be pivoted
+ * (true) or the one with the largest absolute value (false).
+ * @return 0 if no interchange occured, -1 if all are zero (usually signaling
+ * a degeneracy) and positive integer k means that rows ro and k were swapped.
+ */
+int matrix::pivot(unsigned ro, unsigned co, bool symbolic)
+{
+ unsigned k = ro;
+ if (symbolic) {
+ // search first non-zero element in column co beginning at row ro
+ while ((k<row) && (this->m[k*col+co].expand().is_zero()))
+ ++k;
+ } else {
+ // search largest element in column co beginning at row ro
+ GINAC_ASSERT(is_exactly_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_exactly_a<numeric>(this->m[kmax*col+co]));
+ numeric tmp = ex_to<numeric>(this->m[kmax*col+co]);
+ if (abs(tmp) > mmax) {
+ mmax = tmp;
+ k = kmax;
+ }
+ ++kmax;
+ }
+ if (!mmax.is_zero())
+ k = kmax;
+ }
+ if (k==row)
+ // all elements in column co below row ro vanish
+ return -1;
+ if (k==ro)
+ // matrix needs no pivoting
+ return 0;
+ // matrix needs pivoting, so swap rows k and ro
+ ensure_if_modifiable();
+ for (unsigned c=0; c<col; ++c)
+ this->m[k*col+c].swap(this->m[ro*col+c]);
+
+ return k;
+}
+
+ex lst_to_matrix(const lst & l)
+{
+ lst::const_iterator itr, itc;
+
+ // Find number of rows and columns
+ size_t rows = l.nops(), cols = 0;
+ for (itr = l.begin(); itr != l.end(); ++itr) {
+ if (!is_a<lst>(*itr))
+ throw (std::invalid_argument("lst_to_matrix: argument must be a list of lists"));
+ if (itr->nops() > cols)
+ cols = itr->nops();
+ }
+
+ // Allocate and fill matrix
+ matrix &M = *new matrix(rows, cols);
+ M.setflag(status_flags::dynallocated);
+
+ unsigned i;
+ for (itr = l.begin(), i = 0; itr != l.end(); ++itr, ++i) {
+ unsigned j;
+ for (itc = ex_to<lst>(*itr).begin(), j = 0; itc != ex_to<lst>(*itr).end(); ++itc, ++j)
+ M(i, j) = *itc;
+ }
+
+ return M;
+}
+
+ex diag_matrix(const lst & l)
+{
+ lst::const_iterator it;
+ size_t dim = l.nops();
+
+ // Allocate and fill matrix
+ matrix &M = *new matrix(dim, dim);
+ M.setflag(status_flags::dynallocated);
+
+ unsigned i;
+ for (it = l.begin(), i = 0; it != l.end(); ++it, ++i)
+ M(i, i) = *it;
+
+ return M;
+}
+
+ex unit_matrix(unsigned r, unsigned c)
+{
+ matrix &Id = *new matrix(r, c);
+ Id.setflag(status_flags::dynallocated);
+ for (unsigned i=0; i<r && i<c; i++)
+ Id(i,i) = _ex1;
+
+ return Id;
+}
+
+ex symbolic_matrix(unsigned r, unsigned c, const std::string & base_name, const std::string & tex_base_name)
+{
+ matrix &M = *new matrix(r, c);
+ M.setflag(status_flags::dynallocated | status_flags::evaluated);
+
+ bool long_format = (r > 10 || c > 10);
+ bool single_row = (r == 1 || c == 1);
+
+ for (unsigned i=0; i<r; i++) {
+ for (unsigned j=0; j<c; j++) {
+ std::ostringstream s1, s2;
+ s1 << base_name;
+ s2 << tex_base_name << "_{";
+ if (single_row) {
+ if (c == 1) {
+ s1 << i;
+ s2 << i << '}';
+ } else {
+ s1 << j;
+ s2 << j << '}';
+ }
+ } else {
+ if (long_format) {
+ s1 << '_' << i << '_' << j;
+ s2 << i << ';' << j << "}";
+ } else {
+ s1 << i << j;
+ s2 << i << j << '}';
+ }
+ }
+ M(i, j) = symbol(s1.str(), s2.str());
+ }
+ }
+
+ return M;
+}