+ 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 determinant 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::iterator it = this->m.begin();
+ exvector::iterator tmp_n_it = tmp_n.m.begin();
+ exvector::iterator tmp_d_it = tmp_d.m.begin();
+ for (; it!= this->m.end(); ++it, ++tmp_n_it, ++tmp_d_it) {
+ (*tmp_n_it) = (*it).normal().to_rational(srl);
+ (*tmp_d_it) = (*tmp_n_it).denom();
+ (*tmp_n_it) = (*tmp_n_it).numer();
+ }
+
+ 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:
+ it = this->m.begin();
+ tmp_n_it = tmp_n.m.begin();
+ tmp_d_it = tmp_d.m.begin();
+ for (; it!= this->m.end(); ++it, ++tmp_n_it, ++tmp_d_it)
+ (*it) = ((*tmp_n_it)/(*tmp_d_it)).subs(srl);
+
+ return sign;