- // 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)}.
-
- GINAC_ASSERT(det || row==col);
- ensure_if_modifiable();
- if (rows()==1)
- return 1;
-
- int sign = 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(row,col); // for denominators, if needed
- lst srl; // symbol replacement list
- exvector::iterator it = m.begin();
- exvector::iterator tmp_n_it = tmp_n.m.begin();
- exvector::iterator tmp_d_it = tmp_d.m.begin();
- for (; it!= 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();
- }
-
- for (unsigned r1=0; r1<row-1; ++r1) {
- int indx = tmp_n.pivot(r1);
- if (det && indx==-1)
- return 0; // FIXME: what to do if det is false, some day?
- if (indx>0) {
- sign = -sign;
- // rows r1 and indx were swapped, so pivot matrix tmp_d:
- for (unsigned c=0; c<col; ++c)
- tmp_d.m[row*indx+c].swap(tmp_d.m[row*r1+c]);
- }
- if (r1>0) {
- divisor_n = tmp_n.m[(r1-1)*col+(r1-1)].expand();
- divisor_d = tmp_d.m[(r1-1)*col+(r1-1)].expand();
- // save space by deleting no longer needed elements:
- if (det) {
- for (unsigned c=0; c<col; ++c) {
- tmp_n.m[(r1-1)*col+c] = 0;
- tmp_d.m[(r1-1)*col+c] = 1;
- }
- }
- }
- for (unsigned r2=r1+1; r2<row; ++r2) {
- for (unsigned c=r1+1; c<col; ++c) {
- dividend_n = (tmp_n.m[r1*col+r1]*tmp_n.m[r2*col+c]*
- tmp_d.m[r2*col+r1]*tmp_d.m[r1*col+c]
- -tmp_n.m[r2*col+r1]*tmp_n.m[r1*col+c]*
- tmp_d.m[r1*col+r1]*tmp_d.m[r2*col+c]).expand();
- dividend_d = (tmp_d.m[r2*col+r1]*tmp_d.m[r1*col+c]*
- tmp_d.m[r1*col+r1]*tmp_d.m[r2*col+c]).expand();
- bool check = divide(dividend_n, divisor_n,
- tmp_n.m[r2*col+c],true);
- check &= divide(dividend_d, divisor_d,
- tmp_d.m[r2*col+c],true);
- GINAC_ASSERT(check);
- }
- // fill up left hand side.
- for (unsigned c=0; c<=r1; ++c)
- tmp_n.m[r2*col+c] = _ex0();
- }
- }
- // repopulate *this matrix:
- it = m.begin();
- tmp_n_it = tmp_n.m.begin();
- tmp_d_it = tmp_d.m.begin();
- for (; it!= m.end(); ++it, ++tmp_n_it, ++tmp_d_it)
- (*it) = ((*tmp_n_it)/(*tmp_d_it)).subs(srl);
-
- return sign;
+ // 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
+ exmap 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 c0=0; c0<n && r0<m-1; ++c0) {
+ // When trying to find a pivot, we should try a bit harder than expand().
+ // Searching the first non-zero element in-place here instead of calling
+ // pivot() allows us to do no more substitutions and back-substitutions
+ // than are actually necessary.
+ unsigned indx = r0;
+ while ((indx<m) &&
+ (tmp_n[indx*n+c0].subs(srl, subs_options::no_pattern).expand().is_zero()))
+ ++indx;
+ if (indx==m) {
+ // all elements in column c0 below row r0 vanish
+ sign = 0;
+ if (det)
+ return 0;
+ } else {
+ if (indx>r0) {
+ // Matrix needs pivoting, swap rows r0 and indx of tmp_n and tmp_d.
+ sign = -sign;
+ for (unsigned c=c0; c<n; ++c) {
+ tmp_n.m[n*indx+c].swap(tmp_n.m[n*r0+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=c0+1; c<n; ++c) {
+ dividend_n = (tmp_n.m[r0*n+c0]*tmp_n.m[r2*n+c]*
+ tmp_d.m[r2*n+c0]*tmp_d.m[r0*n+c]
+ -tmp_n.m[r2*n+c0]*tmp_n.m[r0*n+c]*
+ tmp_d.m[r0*n+c0]*tmp_d.m[r2*n+c]).expand();
+ dividend_d = (tmp_d.m[r2*n+c0]*tmp_d.m[r0*n+c]*
+ tmp_d.m[r0*n+c0]*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=r0; c<=c0; ++c)
+ tmp_n.m[r2*n+c] = _ex0;
+ }
+ if (c0<n && r0<m-1) {
+ // compute next iteration's divisor
+ divisor_n = tmp_n.m[r0*n+c0].expand();
+ divisor_d = tmp_d.m[r0*n+c0].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;
+ }
+ }
+ // clear remaining rows
+ for (unsigned r=r0+1; r<m; ++r) {
+ for (unsigned c=0; c<n; ++c)
+ tmp_n.m[r*n+c] = _ex0;
+ }
+
+ // 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;