1 /*
2 * GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program; if not, write to the Free Software
16 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
17 */
18
19 #ifndef _GEOM_SL_DETERMINANT_MINOR_H_
20 #define _GEOM_SL_DETERMINANT_MINOR_H_
21
22 #include <map>
23
24
25 namespace Geom { namespace SL {
26
27 /*
28 * determinant_minor
29 * This routine has been taken from the ginac project
30 * and adapted as needed; comments are the original ones.
31 */
32
33 /** Recursive determinant for small matrices having at least one symbolic
34 * entry. The basic algorithm, known as Laplace-expansion, is enhanced by
35 * some bookkeeping to avoid calculation of the same submatrices ("minors")
36 * more than once. According to W.M.Gentleman and S.C.Johnson this algorithm
37 * is better than elimination schemes for matrices of sparse multivariate
38 * polynomials and also for matrices of dense univariate polynomials if the
39 * matrix' dimesion is larger than 7.
40 *
41 * @return the determinant as a new expression (in expanded form)
42 * @see matrix::determinant() */
43
44 template< typename Coeff >
determinant_minor(Matrix<Coeff> const & M)45 Coeff determinant_minor(Matrix<Coeff> const& M)
46 {
47 assert(M.rows() == M.columns());
48 // for small matrices the algorithm does not make any sense:
49 const unsigned int n = M.columns();
50 if (n == 1)
51 return M(0,0);
52 if (n == 2)
53 return (M(0,0) * M(1,1) - M(0,1) * M(1,0));
54 if (n == 3)
55 return ( M(0,0)*M(1,1)*M(2,2) + M(0,2)*M(1,0)*M(2,1)
56 + M(0,1)*M(1,2)*M(2,0) - M(0,2)*M(1,1)*M(2,0)
57 - M(0,0)*M(1,2)*M(2,1) - M(0,1)*M(1,0)*M(2,2) );
58
59 // This algorithm can best be understood by looking at a naive
60 // implementation of Laplace-expansion, like this one:
61 // ex det;
62 // matrix minorM(this->rows()-1,this->cols()-1);
63 // for (unsigned r1=0; r1<this->rows(); ++r1) {
64 // // shortcut if element(r1,0) vanishes
65 // if (m[r1*col].is_zero())
66 // continue;
67 // // assemble the minor matrix
68 // for (unsigned r=0; r<minorM.rows(); ++r) {
69 // for (unsigned c=0; c<minorM.cols(); ++c) {
70 // if (r<r1)
71 // minorM(r,c) = m[r*col+c+1];
72 // else
73 // minorM(r,c) = m[(r+1)*col+c+1];
74 // }
75 // }
76 // // recurse down and care for sign:
77 // if (r1%2)
78 // det -= m[r1*col] * minorM.determinant_minor();
79 // else
80 // det += m[r1*col] * minorM.determinant_minor();
81 // }
82 // return det.expand();
83 // What happens is that while proceeding down many of the minors are
84 // computed more than once. In particular, there are binomial(n,k)
85 // kxk minors and each one is computed factorial(n-k) times. Therefore
86 // it is reasonable to store the results of the minors. We proceed from
87 // right to left. At each column c we only need to retrieve the minors
88 // calculated in step c-1. We therefore only have to store at most
89 // 2*binomial(n,n/2) minors.
90
91 // Unique flipper counter for partitioning into minors
92 std::vector<unsigned int> Pkey;
93 Pkey.reserve(n);
94 // key for minor determinant (a subpartition of Pkey)
95 std::vector<unsigned int> Mkey;
96 Mkey.reserve(n-1);
97 // we store our subminors in maps, keys being the rows they arise from
98 typedef typename std::map<std::vector<unsigned>, Coeff> Rmap;
99 typedef typename std::map<std::vector<unsigned>, Coeff>::value_type Rmap_value;
100 Rmap A;
101 Rmap B;
102 Coeff det;
103 // initialize A with last column:
104 for (unsigned int r = 0; r < n; ++r)
105 {
106 Pkey.erase(Pkey.begin(),Pkey.end());
107 Pkey.push_back(r);
108 A.insert(Rmap_value(Pkey,M(r,n-1)));
109 }
110 // proceed from right to left through matrix
111 for (int c = n-2; c >= 0; --c)
112 {
113 Pkey.erase(Pkey.begin(),Pkey.end()); // don't change capacity
114 Mkey.erase(Mkey.begin(),Mkey.end());
115 for (unsigned int i = 0; i < n-c; ++i)
116 Pkey.push_back(i);
117 unsigned int fc = 0; // controls logic for our strange flipper counter
118 do
119 {
120 det = Geom::SL::zero<Coeff>()();
121 for (unsigned int r = 0; r < n-c; ++r)
122 {
123 // maybe there is nothing to do?
124 if (M(Pkey[r], c).is_zero())
125 continue;
126 // create the sorted key for all possible minors
127 Mkey.erase(Mkey.begin(),Mkey.end());
128 for (unsigned int i = 0; i < n-c; ++i)
129 if (i != r)
130 Mkey.push_back(Pkey[i]);
131 // Fetch the minors and compute the new determinant
132 if (r % 2)
133 det -= M(Pkey[r],c)*A[Mkey];
134 else
135 det += M(Pkey[r],c)*A[Mkey];
136 }
137 // store the new determinant at its place in B:
138 if (!det.is_zero())
139 B.insert(Rmap_value(Pkey,det));
140 // increment our strange flipper counter
141 for (fc = n-c; fc > 0; --fc)
142 {
143 ++Pkey[fc-1];
144 if (Pkey[fc-1]<fc+c)
145 break;
146 }
147 if (fc < n-c && fc > 0)
148 for (unsigned int j = fc; j < n-c; ++j)
149 Pkey[j] = Pkey[j-1]+1;
150 } while(fc);
151 // next column, so change the role of A and B:
152 A.swap(B);
153 B.clear();
154 }
155
156 return det;
157 }
158
159
160
161 } /*end namespace Geom*/ } /*end namespace SL*/
162
163 #endif // _GEOM_SL_DETERMINANT_MINOR_H_
164
165
166 /*
167 Local Variables:
168 mode:c++
169 c-file-style:"stroustrup"
170 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
171 indent-tabs-mode:nil
172 fill-column:99
173 End:
174 */
175 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
176