[bugfix] parser::parse_literal_expr(): don't forget to consume the token...
[ginac.git] / check / check_matrices.cpp
1 /** @file check_matrices.cpp
2  *
3  *  Here we test manipulations on GiNaC's symbolic matrices.  They are a
4  *  well-tried resource for cross-checking the underlying symbolic
5  *  manipulations. */
6
7 /*
8  *  GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
9  *
10  *  This program is free software; you can redistribute it and/or modify
11  *  it under the terms of the GNU General Public License as published by
12  *  the Free Software Foundation; either version 2 of the License, or
13  *  (at your option) any later version.
14  *
15  *  This program is distributed in the hope that it will be useful,
16  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18  *  GNU General Public License for more details.
19  *
20  *  You should have received a copy of the GNU General Public License
21  *  along with this program; if not, write to the Free Software
22  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
23  */
24
25 #include <iostream>
26 #include <cstdlib> // rand(), RAND_MAX
27 #include "ginac.h"
28 using namespace std;
29 using namespace GiNaC;
30
31 extern const ex 
32 sparse_tree(const symbol & x, const symbol & y, const symbol & z,
33             int level, bool trig = false, bool rational = true,
34             bool complex = false);
35 extern const ex 
36 dense_univariate_poly(const symbol & x, unsigned degree);
37
38 /* determinants of some sparse symbolic matrices with coefficients in
39  * an integral domain. */
40 static unsigned integdom_matrix_determinants()
41 {
42         unsigned result = 0;
43         symbol a("a");
44         
45         for (unsigned size=3; size<22; ++size) {
46                 matrix A(size,size);
47                 // populate one element in each row:
48                 for (unsigned r=0; r<size-1; ++r)
49                         A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
50                 // set the last row to a linear combination of two other lines
51                 // to guarantee that the determinant is zero:
52                 for (unsigned c=0; c<size; ++c)
53                         A.set(size-1,c,A(0,c)-A(size-2,c));
54                 if (!A.determinant().is_zero()) {
55                         clog << "Determinant of " << size << "x" << size << " matrix "
56                              << endl << A << endl
57                              << "was not found to vanish!" << endl;
58                         ++result;
59                 }
60         }
61         
62         return result;
63 }
64
65 /* determinants of some symbolic matrices with multivariate rational function
66  * coefficients. */
67 static unsigned rational_matrix_determinants()
68 {
69         unsigned result = 0;
70         symbol a("a"), b("b"), c("c");
71         
72         for (unsigned size=3; size<9; ++size) {
73                 matrix A(size,size);
74                 for (unsigned r=0; r<size-1; ++r) {
75                         // populate one or two elements in each row:
76                         for (unsigned ec=0; ec<2; ++ec) {
77                                 ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
78                                 ex denom;
79                                 do {
80                                         denom = sparse_tree(a, b, c, rand()%2, false, false, false);
81                                 } while (denom.is_zero());
82                                 A.set(r,unsigned(rand()%size),numer/denom);
83                         }
84                 }
85                 // set the last row to a linear combination of two other lines
86                 // to guarantee that the determinant is zero:
87                 for (unsigned co=0; co<size; ++co)
88                         A.set(size-1,co,A(0,co)-A(size-2,co));
89                 if (!A.determinant().is_zero()) {
90                         clog << "Determinant of " << size << "x" << size << " matrix "
91                              << endl << A << endl
92                              << "was not found to vanish!" << endl;
93                         ++result;
94                 }
95         }
96         
97         return result;
98 }
99
100 /* Some quite funny determinants with functions and stuff like that inside. */
101 static unsigned funny_matrix_determinants()
102 {
103         unsigned result = 0;
104         symbol a("a"), b("b"), c("c");
105         
106         for (unsigned size=3; size<8; ++size) {
107                 matrix A(size,size);
108                 for (unsigned co=0; co<size-1; ++co) {
109                         // populate one or two elements in each row:
110                         for (unsigned ec=0; ec<2; ++ec) {
111                                 ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
112                                 ex denom;
113                                 do {
114                                         denom = sparse_tree(a, b, c, rand()%2, false, true, false);
115                                 } while (denom.is_zero());
116                                 A.set(unsigned(rand()%size),co,numer/denom);
117                         }
118                 }
119                 // set the last column to a linear combination of two other columns
120                 // to guarantee that the determinant is zero:
121                 for (unsigned ro=0; ro<size; ++ro)
122                         A.set(ro,size-1,A(ro,0)-A(ro,size-2));
123                 if (!A.determinant().is_zero()) {
124                         clog << "Determinant of " << size << "x" << size << " matrix "
125                              << endl << A << endl
126                              << "was not found to vanish!" << endl;
127                         ++result;
128                 }
129         }
130         
131         return result;
132 }
133
134 /* compare results from different determinant algorithms.*/
135 static unsigned compare_matrix_determinants()
136 {
137         unsigned result = 0;
138         symbol a("a");
139         
140         for (unsigned size=2; size<8; ++size) {
141                 matrix A(size,size);
142                 for (unsigned co=0; co<size; ++co) {
143                         for (unsigned ro=0; ro<size; ++ro) {
144                                 // populate some elements
145                                 ex elem = 0;
146                                 if (rand()%(size/2) == 0)
147                                         elem = sparse_tree(a, a, a, rand()%3, false, true, false);
148                                 A.set(ro,co,elem);
149                         }
150                 }
151                 ex det_gauss = A.determinant(determinant_algo::gauss);
152                 ex det_laplace = A.determinant(determinant_algo::laplace);
153                 ex det_divfree = A.determinant(determinant_algo::divfree);
154                 ex det_bareiss = A.determinant(determinant_algo::bareiss);
155                 if ((det_gauss-det_laplace).normal() != 0 ||
156                         (det_bareiss-det_laplace).normal() != 0 ||
157                         (det_divfree-det_laplace).normal() != 0) {
158                         clog << "Determinant of " << size << "x" << size << " matrix "
159                              << endl << A << endl
160                              << "is inconsistent between different algorithms:" << endl
161                              << "Gauss elimination:   " << det_gauss << endl
162                              << "Minor elimination:   " << det_laplace << endl
163                              << "Division-free elim.: " << det_divfree << endl
164                              << "Fraction-free elim.: " << det_bareiss << endl;
165                         ++result;
166                 }
167         }
168         
169         return result;
170 }
171
172 static unsigned symbolic_matrix_inverse()
173 {
174         unsigned result = 0;
175         symbol a("a"), b("b"), c("c");
176         
177         for (unsigned size=2; size<6; ++size) {
178                 matrix A(size,size);
179                 do {
180                         for (unsigned co=0; co<size; ++co) {
181                                 for (unsigned ro=0; ro<size; ++ro) {
182                                         // populate some elements
183                                         ex elem = 0;
184                                         if (rand()%(size/2) == 0)
185                                                 elem = sparse_tree(a, b, c, rand()%2, false, true, false);
186                                         A.set(ro,co,elem);
187                                 }
188                         }
189                 } while (A.determinant() == 0);
190                 matrix B = A.inverse();
191                 matrix C = A.mul(B);
192                 bool ok = true;
193                 for (unsigned ro=0; ro<size; ++ro)
194                         for (unsigned co=0; co<size; ++co)
195                                 if (C(ro,co).normal() != (ro==co?1:0))
196                                         ok = false;
197                 if (!ok) {
198                         clog << "Inverse of " << size << "x" << size << " matrix "
199                              << endl << A << endl
200                              << "erroneously returned: "
201                              << endl << B << endl;
202                         ++result;
203                 }
204         }
205         
206         return result;
207 }
208
209 unsigned check_matrices()
210 {
211         unsigned result = 0;
212         
213         cout << "checking symbolic matrix manipulations" << flush;
214         
215         result += integdom_matrix_determinants();  cout << '.' << flush;
216         result += rational_matrix_determinants();  cout << '.' << flush;
217         result += funny_matrix_determinants();  cout << '.' << flush;
218         result += compare_matrix_determinants();  cout << '.' << flush;
219         result += symbolic_matrix_inverse();  cout << '.' << flush;
220         
221         return result;
222 }
223
224 int main(int argc, char** argv)
225 {
226         return check_matrices();
227 }