+ // 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_symbolic_minor();
+ // else
+ // det += m[r1*col] * minorM.determinant_symbolic_minor();
+ // }
+ // return det;
+ // 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.
+
+ // 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, B;
+ ex det;
+ vector<unsigned> Pkey; // Unique flipper counter for partitioning into minors
+ Pkey.reserve(this->col);
+ vector<unsigned> Mkey; // key for minor determinant (a subpartition of Pkey)
+ Mkey.reserve(this->col-1);
+ // 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 {
+ A.insert(Rmap_value(Pkey,_ex0()));
+ 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];