- // 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();
-
- // This algorithm can best be understood by looking at a naive
- // implementation of Laplace-expansion, like this one:
- // 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();
- // else
- // det += m[r1*col] * minorM.determinant_minor();
- // }
- // return det.expand();
- // What happens is that while proceeding down many of the minors are
- // computed more than once. In particular, there are binomial(n,k)
- // kxk minors and each one is computed factorial(n-k) times. Therefore
- // it is reasonable to store the results of the minors. We proceed from
- // right to left. At each column c we only need to retrieve the minors
- // calculated in step c-1. We therefore only have to store at most
- // 2*binomial(n,n/2) minors.
-
- // Unique flipper counter for partitioning into minors
- vector<unsigned> Pkey;
- Pkey.reserve(this->col);
- // key for minor determinant (a subpartition of Pkey)
- vector<unsigned> Mkey;
- Mkey.reserve(this->col-1);
- // we store our subminors in maps, keys being the rows they arise from
- typedef map<vector<unsigned>,class ex> Rmap;
- typedef map<vector<unsigned>,class ex>::value_type Rmap_value;
- Rmap A;
- Rmap B;
- ex det;
- // initialize A with last column:
- for (unsigned r=0; r<this->col; ++r) {
- Pkey.erase(Pkey.begin(),Pkey.end());
- Pkey.push_back(r);
- A.insert(Rmap_value(Pkey,m[this->col*r+this->col-1]));
- }
- // proceed from right to left through matrix
- for (int c=this->col-2; c>=0; --c) {
- Pkey.erase(Pkey.begin(),Pkey.end()); // don't change capacity
- Mkey.erase(Mkey.begin(),Mkey.end());
- for (unsigned i=0; i<this->col-c; ++i)
- Pkey.push_back(i);
- unsigned fc = 0; // controls logic for our strange flipper counter
- do {
- det = _ex0();
- for (unsigned r=0; r<this->col-c; ++r) {
- // maybe there is nothing to do?
- if (m[Pkey[r]*this->col+c].is_zero())
- continue;
- // create the sorted key for all possible minors
- Mkey.erase(Mkey.begin(),Mkey.end());
- for (unsigned i=0; i<this->col-c; ++i)
- if (i!=r)
- Mkey.push_back(Pkey[i]);
- // Fetch the minors and compute the new determinant
- if (r%2)
- det -= m[Pkey[r]*this->col+c]*A[Mkey];
- else
- det += m[Pkey[r]*this->col+c]*A[Mkey];
- }
- // prevent build-up of deep nesting of expressions saves time:
- det = det.expand();
- // store the new determinant at its place in B:
- if (!det.is_zero())
- B.insert(Rmap_value(Pkey,det));
- // increment our strange flipper counter
- for (fc=this->col-c; fc>0; --fc) {
- ++Pkey[fc-1];
- if (Pkey[fc-1]<fc+c)
- break;
- }
- if (fc<this->col-c)
- for (unsigned j=fc; j<this->col-c; ++j)
- Pkey[j] = Pkey[j-1]+1;
- } while(fc);
- // next column, so change the role of A and B:
- A = B;
- B.clear();
- }
-
- return det;
-}
-
-/** 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)
-{
- 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 computed by using fraction free elimination. This
- * routine is only called internally by matrix::determinant().
- *
- * @param normalize may be set to false only in integral domains. */
-ex matrix::determinant_bareiss(bool normalize) const
-{
- if (rows()==1)
- return m[0];
-
- int sign = 1;
- ex divisor = 1;
- ex dividend;
-
- // 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);
- }
-
- 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];
+ // for small matrices the algorithm does not make any sense:
+ const unsigned n = this->cols();
+ if (n==1)
+ return m[0].expand();
+ if (n==2)
+ return (m[0]*m[3]-m[2]*m[1]).expand();
+ if (n==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();
+
+ // This algorithm can best be understood by looking at a naive
+ // implementation of Laplace-expansion, like this one:
+ // ex det;
+ // matrix minorM(this->rows()-1,this->cols()-1);
+ // for (unsigned r1=0; r1<this->rows(); ++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(r,c) = m[r*col+c+1];
+ // else
+ // minorM(r,c) = m[(r+1)*col+c+1];
+ // }
+ // }
+ // // recurse down and care for sign:
+ // if (r1%2)
+ // det -= m[r1*col] * minorM.determinant_minor();
+ // else
+ // det += m[r1*col] * minorM.determinant_minor();
+ // }
+ // return det.expand();
+ // What happens is that while proceeding down many of the minors are
+ // computed more than once. In particular, there are binomial(n,k)
+ // kxk minors and each one is computed factorial(n-k) times. Therefore
+ // it is reasonable to store the results of the minors. We proceed from
+ // right to left. At each column c we only need to retrieve the minors
+ // calculated in step c-1. We therefore only have to store at most
+ // 2*binomial(n,n/2) minors.
+
+ // Unique flipper counter for partitioning into minors
+ std::vector<unsigned> Pkey;
+ Pkey.reserve(n);
+ // key for minor determinant (a subpartition of Pkey)
+ std::vector<unsigned> Mkey;
+ Mkey.reserve(n-1);
+ // we store our subminors in maps, keys being the rows they arise from
+ typedef std::map<std::vector<unsigned>,class ex> Rmap;
+ typedef std::map<std::vector<unsigned>,class ex>::value_type Rmap_value;
+ Rmap A;
+ Rmap B;
+ ex det;
+ // initialize A with last column:
+ for (unsigned r=0; r<n; ++r) {
+ Pkey.erase(Pkey.begin(),Pkey.end());
+ Pkey.push_back(r);
+ A.insert(Rmap_value(Pkey,m[n*(r+1)-1]));
+ }
+ // proceed from right to left through matrix
+ for (int c=n-2; c>=0; --c) {
+ Pkey.erase(Pkey.begin(),Pkey.end()); // don't change capacity
+ Mkey.erase(Mkey.begin(),Mkey.end());
+ for (unsigned i=0; i<n-c; ++i)
+ Pkey.push_back(i);
+ unsigned fc = 0; // controls logic for our strange flipper counter
+ do {
+ det = _ex0;
+ for (unsigned r=0; r<n-c; ++r) {
+ // maybe there is nothing to do?
+ if (m[Pkey[r]*n+c].is_zero())
+ continue;
+ // create the sorted key for all possible minors
+ Mkey.erase(Mkey.begin(),Mkey.end());
+ for (unsigned i=0; i<n-c; ++i)
+ if (i!=r)
+ Mkey.push_back(Pkey[i]);
+ // Fetch the minors and compute the new determinant
+ if (r%2)
+ det -= m[Pkey[r]*n+c]*A[Mkey];
+ else
+ det += m[Pkey[r]*n+c]*A[Mkey];
+ }
+ // prevent build-up of deep nesting of expressions saves time:
+ det = det.expand();
+ // store the new determinant at its place in B:
+ if (!det.is_zero())
+ B.insert(Rmap_value(Pkey,det));
+ // increment our strange flipper counter
+ for (fc=n-c; fc>0; --fc) {
+ ++Pkey[fc-1];
+ if (Pkey[fc-1]<fc+c)
+ break;
+ }
+ if (fc<n-c && fc>0)
+ for (unsigned j=fc; j<n-c; ++j)
+ Pkey[j] = Pkey[j-1]+1;
+ } while(fc);
+ // next column, so change the role of A and B:
+ A = B;
+ B.clear();
+ }
+
+ return det;