-#endif // def DOASSERT
-
- // assemble solution
- matrix sol(n,1);
- int last_assigned_sol=n+1;
- for (int r=m; r>0; --r) {
- int first_non_zero=1;
- while ((first_non_zero<=n)&&(a.ffe_get(r,first_non_zero).is_zero())) {
- first_non_zero++;
- }
- if (first_non_zero>n) {
- // row consists only of zeroes, corresponding rhs must be 0 as well
- if (!b.ffe_get(r,1).is_zero()) {
- throw (std::runtime_error("matrix::fraction_free_elim(): singular matrix"));
- }
- } else {
- // assign solutions for vars between first_non_zero+1 and
- // last_assigned_sol-1: free parameters
- for (int c=first_non_zero+1; c<=last_assigned_sol-1; ++c) {
- sol.ffe_set(c,1,vars.ffe_get(c,1));
+
+ return sol;
+}
+
+
+// protected
+
+/** 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")
+ * more than once. According to W.M.Gentleman and S.C.Johnson this algorithm
+ * is better than elimination schemes for matrices of sparse multivariate
+ * polynomials and also for matrices of dense univariate polynomials if the
+ * matrix' dimesion is larger than 7.
+ *
+ * @return the determinant as a new expression (in expanded form)
+ * @see matrix::determinant() */
+ex matrix::determinant_minor(void) const
+{
+ // 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.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
+ 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];