Happy New Year!
[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-2019 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 "ginac.h"
26 using namespace GiNaC;
27
28 #include <cstdlib> // for rand(), RAND_MAX
29 #include <iostream>
30 using namespace std;
31
32 extern const ex 
33 sparse_tree(const symbol & x, const symbol & y, const symbol & z,
34             int level, bool trig = false, bool rational = true,
35             bool complex = false);
36 extern const ex 
37 dense_univariate_poly(const symbol & x, unsigned degree);
38
39 /* determinants of some sparse symbolic matrices with coefficients in
40  * an integral domain. */
41 static unsigned integdom_matrix_determinants()
42 {
43         unsigned result = 0;
44         symbol a("a");
45         
46         for (unsigned size=3; size<22; ++size) {
47                 matrix A(size,size);
48                 // populate one element in each row:
49                 for (unsigned r=0; r<size-1; ++r)
50                         A.set(r,unsigned(rand()%size),dense_univariate_poly(a,5));
51                 // set the last row to a linear combination of two other lines
52                 // to guarantee that the determinant is zero:
53                 for (unsigned c=0; c<size; ++c)
54                         A.set(size-1,c,A(0,c)-A(size-2,c));
55                 if (!A.determinant().is_zero()) {
56                         clog << "Determinant of " << size << "x" << size << " matrix "
57                              << endl << A << endl
58                              << "was not found to vanish!" << endl;
59                         ++result;
60                 }
61         }
62         
63         return result;
64 }
65
66 /* determinants of some symbolic matrices with multivariate rational function
67  * coefficients. */
68 static unsigned rational_matrix_determinants()
69 {
70         unsigned result = 0;
71         symbol a("a"), b("b"), c("c");
72         
73         for (unsigned size=3; size<9; ++size) {
74                 matrix A(size,size);
75                 for (unsigned r=0; r<size-1; ++r) {
76                         // populate one or two elements in each row:
77                         for (unsigned ec=0; ec<2; ++ec) {
78                                 ex numer = sparse_tree(a, b, c, 1+rand()%4, false, false, false);
79                                 ex denom;
80                                 do {
81                                         denom = sparse_tree(a, b, c, rand()%2, false, false, false);
82                                 } while (denom.is_zero());
83                                 A.set(r,unsigned(rand()%size),numer/denom);
84                         }
85                 }
86                 // set the last row to a linear combination of two other lines
87                 // to guarantee that the determinant is zero:
88                 for (unsigned co=0; co<size; ++co)
89                         A.set(size-1,co,A(0,co)-A(size-2,co));
90                 if (!A.determinant().is_zero()) {
91                         clog << "Determinant of " << size << "x" << size << " matrix "
92                              << endl << A << endl
93                              << "was not found to vanish!" << endl;
94                         ++result;
95                 }
96         }
97         
98         return result;
99 }
100
101 /* Some quite funny determinants with functions and stuff like that inside. */
102 static unsigned funny_matrix_determinants()
103 {
104         unsigned result = 0;
105         symbol a("a"), b("b"), c("c");
106         
107         for (unsigned size=3; size<8; ++size) {
108                 matrix A(size,size);
109                 for (unsigned co=0; co<size-1; ++co) {
110                         // populate one or two elements in each row:
111                         for (unsigned ec=0; ec<2; ++ec) {
112                                 ex numer = sparse_tree(a, b, c, 1+rand()%3, true, true, false);
113                                 ex denom;
114                                 do {
115                                         denom = sparse_tree(a, b, c, rand()%2, false, true, false);
116                                 } while (denom.is_zero());
117                                 A.set(unsigned(rand()%size),co,numer/denom);
118                         }
119                 }
120                 // set the last column to a linear combination of two other columns
121                 // to guarantee that the determinant is zero:
122                 for (unsigned ro=0; ro<size; ++ro)
123                         A.set(ro,size-1,A(ro,0)-A(ro,size-2));
124                 if (!A.determinant().is_zero()) {
125                         clog << "Determinant of " << size << "x" << size << " matrix "
126                              << endl << A << endl
127                              << "was not found to vanish!" << endl;
128                         ++result;
129                 }
130         }
131         
132         return result;
133 }
134
135 /* compare results from different determinant algorithms.*/
136 static unsigned compare_matrix_determinants()
137 {
138         unsigned result = 0;
139         symbol a("a");
140         
141         for (unsigned size=2; size<8; ++size) {
142                 matrix A(size,size);
143                 for (unsigned co=0; co<size; ++co) {
144                         for (unsigned ro=0; ro<size; ++ro) {
145                                 // populate some elements
146                                 ex elem = 0;
147                                 if (rand()%(size/2) == 0)
148                                         elem = sparse_tree(a, a, a, rand()%3, false, true, false);
149                                 A.set(ro,co,elem);
150                         }
151                 }
152                 ex det_gauss = A.determinant(determinant_algo::gauss);
153                 ex det_laplace = A.determinant(determinant_algo::laplace);
154                 ex det_divfree = A.determinant(determinant_algo::divfree);
155                 ex det_bareiss = A.determinant(determinant_algo::bareiss);
156                 if ((det_gauss-det_laplace).normal() != 0 ||
157                         (det_bareiss-det_laplace).normal() != 0 ||
158                         (det_divfree-det_laplace).normal() != 0) {
159                         clog << "Determinant of " << size << "x" << size << " matrix "
160                              << endl << A << endl
161                              << "is inconsistent between different algorithms:" << endl
162                              << "Gauss elimination:   " << det_gauss << endl
163                              << "Minor elimination:   " << det_laplace << endl
164                              << "Division-free elim.: " << det_divfree << endl
165                              << "Fraction-free elim.: " << det_bareiss << endl;
166                         ++result;
167                 }
168         }
169         
170         return result;
171 }
172
173 static unsigned symbolic_matrix_inverse()
174 {
175         unsigned result = 0;
176         symbol a("a"), b("b"), c("c");
177         
178         for (unsigned size=2; size<6; ++size) {
179                 matrix A(size,size);
180                 do {
181                         for (unsigned co=0; co<size; ++co) {
182                                 for (unsigned ro=0; ro<size; ++ro) {
183                                         // populate some elements
184                                         ex elem = 0;
185                                         if (rand()%(size/2) == 0)
186                                                 elem = sparse_tree(a, b, c, rand()%2, false, true, false);
187                                         A.set(ro,co,elem);
188                                 }
189                         }
190                 } while (A.determinant() == 0);
191                 matrix B = A.inverse();
192                 matrix C = A.mul(B);
193                 bool ok = true;
194                 for (unsigned ro=0; ro<size; ++ro)
195                         for (unsigned co=0; co<size; ++co)
196                                 if (C(ro,co).normal() != (ro==co?1:0))
197                                         ok = false;
198                 if (!ok) {
199                         clog << "Inverse of " << size << "x" << size << " matrix "
200                              << endl << A << endl
201                              << "erroneously returned: "
202                              << endl << B << endl;
203                         ++result;
204                 }
205         }
206         
207         return result;
208 }
209
210 unsigned check_matrices()
211 {
212         unsigned result = 0;
213         
214         cout << "checking symbolic matrix manipulations" << flush;
215         
216         result += integdom_matrix_determinants();  cout << '.' << flush;
217         result += rational_matrix_determinants();  cout << '.' << flush;
218         result += funny_matrix_determinants();  cout << '.' << flush;
219         result += compare_matrix_determinants();  cout << '.' << flush;
220         result += symbolic_matrix_inverse();  cout << '.' << flush;
221         
222         return result;
223 }
224
225 int main(int argc, char** argv)
226 {
227         return check_matrices();
228 }