G_numeric: put convergence/acceleration transofrmations into helper functions.
[ginac.git] / check / exam_matrices.cpp
1 /** @file exam_matrices.cpp
2  *
3  *  Here we examine manipulations on GiNaC's symbolic matrices. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
7  *
8  *  This program is free software; you can redistribute it and/or modify
9  *  it under the terms of the GNU General Public License as published by
10  *  the Free Software Foundation; either version 2 of the License, or
11  *  (at your option) any later version.
12  *
13  *  This program is distributed in the hope that it will be useful,
14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  *  GNU General Public License for more details.
17  *
18  *  You should have received a copy of the GNU General Public License
19  *  along with this program; if not, write to the Free Software
20  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
21  */
22
23 #include <stdexcept>
24 #include <iostream>
25 #include "ginac.h"
26 using namespace std;
27 using namespace GiNaC;
28
29 static unsigned matrix_determinants()
30 {
31         unsigned result = 0;
32         ex det;
33         matrix m1(1,1), m2(2,2), m3(3,3), m4(4,4);
34         symbol a("a"), b("b"), c("c");
35         symbol d("d"), e("e"), f("f");
36         symbol g("g"), h("h"), i("i");
37         
38         // check symbolic trivial matrix determinant
39         m1.set(0,0,a);
40         det = m1.determinant();
41         if (det != a) {
42                 clog << "determinant of 1x1 matrix " << m1
43                      << " erroneously returned " << det << endl;
44                 ++result;
45         }
46         
47         // check generic dense symbolic 2x2 matrix determinant
48         m2.set(0,0,a).set(0,1,b);
49         m2.set(1,0,c).set(1,1,d);
50         det = m2.determinant();
51         if (det != (a*d-b*c)) {
52                 clog << "determinant of 2x2 matrix " << m2
53                      << " erroneously returned " << det << endl;
54                 ++result;
55         }
56         
57         // check generic dense symbolic 3x3 matrix determinant
58         m3.set(0,0,a).set(0,1,b).set(0,2,c);
59         m3.set(1,0,d).set(1,1,e).set(1,2,f);
60         m3.set(2,0,g).set(2,1,h).set(2,2,i);
61         det = m3.determinant();
62         if (det != (a*e*i - a*f*h - d*b*i + d*c*h + g*b*f - g*c*e)) {
63                 clog << "determinant of 3x3 matrix " << m3
64                      << " erroneously returned " << det << endl;
65                 ++result;
66         }
67         
68         // check dense numeric 3x3 matrix determinant
69         m3.set(0,0,numeric(0)).set(0,1,numeric(-1)).set(0,2,numeric(3));
70         m3.set(1,0,numeric(3)).set(1,1,numeric(-2)).set(1,2,numeric(2));
71         m3.set(2,0,numeric(3)).set(2,1,numeric(4)).set(2,2,numeric(-2));
72         det = m3.determinant();
73         if (det != 42) {
74                 clog << "determinant of 3x3 matrix " << m3
75                      << " erroneously returned " << det << endl;
76                 ++result;
77         }
78         
79         // check dense symbolic 2x2 matrix determinant
80         m2.set(0,0,a/(a-b)).set(0,1,1);
81         m2.set(1,0,b/(a-b)).set(1,1,1);
82         det = m2.determinant();
83         if (det != 1) {
84                 if (det.normal() == 1)  // only half wrong
85                         clog << "determinant of 2x2 matrix " << m2
86                              << " was returned unnormalized as " << det << endl;
87                 else  // totally wrong
88                         clog << "determinant of 2x2 matrix " << m2
89                              << " erroneously returned " << det << endl;
90                 ++result;
91         }
92         
93         // check sparse symbolic 4x4 matrix determinant
94         m4.set(0,1,a).set(1,0,b).set(3,2,c).set(2,3,d);
95         det = m4.determinant();
96         if (det != a*b*c*d) {
97                 clog << "determinant of 4x4 matrix " << m4
98                      << " erroneously returned " << det << endl;
99                 ++result;
100         }
101         
102         // check characteristic polynomial
103         m3.set(0,0,a).set(0,1,-2).set(0,2,2);
104         m3.set(1,0,3).set(1,1,a-1).set(1,2,2);
105         m3.set(2,0,3).set(2,1,4).set(2,2,a-3);
106         ex p = m3.charpoly(a);
107         if (p != 0) {
108                 clog << "charpoly of 3x3 matrix " << m3
109                      << " erroneously returned " << p << endl;
110                 ++result;
111         }
112         
113         return result;
114 }
115
116 static unsigned matrix_invert1()
117 {
118         unsigned result = 0;
119         matrix m(1,1);
120         symbol a("a");
121         
122         m.set(0,0,a);
123         matrix m_i = m.inverse();
124         
125         if (m_i(0,0) != pow(a,-1)) {
126                 clog << "inversion of 1x1 matrix " << m
127                      << " erroneously returned " << m_i << endl;
128                 ++result;
129         }
130         
131         return result;
132 }
133
134 static unsigned matrix_invert2()
135 {
136         unsigned result = 0;
137         matrix m(2,2);
138         symbol a("a"), b("b"), c("c"), d("d");
139         m.set(0,0,a).set(0,1,b);
140         m.set(1,0,c).set(1,1,d);
141         matrix m_i = m.inverse();
142         ex det = m.determinant();
143         
144         if ((normal(m_i(0,0)*det) != d) ||
145                 (normal(m_i(0,1)*det) != -b) ||
146                 (normal(m_i(1,0)*det) != -c) ||
147                 (normal(m_i(1,1)*det) != a)) {
148                 clog << "inversion of 2x2 matrix " << m
149                      << " erroneously returned " << m_i << endl;
150                 ++result;
151         }
152         
153         return result;
154 }
155
156 static unsigned matrix_invert3()
157 {
158         unsigned result = 0;
159         matrix m(3,3);
160         symbol a("a"), b("b"), c("c");
161         symbol d("d"), e("e"), f("f");
162         symbol g("g"), h("h"), i("i");
163         m.set(0,0,a).set(0,1,b).set(0,2,c);
164         m.set(1,0,d).set(1,1,e).set(1,2,f);
165         m.set(2,0,g).set(2,1,h).set(2,2,i);
166         matrix m_i = m.inverse();
167         ex det = m.determinant();
168         
169         if ((normal(m_i(0,0)*det) != (e*i-f*h)) ||
170             (normal(m_i(0,1)*det) != (c*h-b*i)) ||
171             (normal(m_i(0,2)*det) != (b*f-c*e)) ||
172             (normal(m_i(1,0)*det) != (f*g-d*i)) ||
173             (normal(m_i(1,1)*det) != (a*i-c*g)) ||
174             (normal(m_i(1,2)*det) != (c*d-a*f)) ||
175             (normal(m_i(2,0)*det) != (d*h-e*g)) ||
176             (normal(m_i(2,1)*det) != (b*g-a*h)) ||
177             (normal(m_i(2,2)*det) != (a*e-b*d))) {
178                 clog << "inversion of 3x3 matrix " << m
179                      << " erroneously returned " << m_i << endl;
180                 ++result;
181         }
182         
183         return result;
184 }
185
186 static unsigned matrix_solve2()
187 {
188         // check the solution of the multiple system A*X = B:
189         //       [ 1  2 -1 ] [ x0 y0 ]   [ 4 0 ]
190         //       [ 1  4 -2 ]*[ x1 y1 ] = [ 7 0 ]
191         //       [ a -2  2 ] [ x2 y2 ]   [ a 4 ]
192         unsigned result = 0;
193         symbol a("a");
194         symbol x0("x0"), x1("x1"), x2("x2");
195         symbol y0("y0"), y1("y1"), y2("y2");
196         matrix A(3,3);
197         A.set(0,0,1).set(0,1,2).set(0,2,-1);
198         A.set(1,0,1).set(1,1,4).set(1,2,-2);
199         A.set(2,0,a).set(2,1,-2).set(2,2,2);
200         matrix B(3,2);
201         B.set(0,0,4).set(1,0,7).set(2,0,a);
202         B.set(0,1,0).set(1,1,0).set(2,1,4);
203         matrix X(3,2);
204         X.set(0,0,x0).set(1,0,x1).set(2,0,x2);
205         X.set(0,1,y0).set(1,1,y1).set(2,1,y2);
206         matrix cmp(3,2);
207         cmp.set(0,0,1).set(1,0,3).set(2,0,3);
208         cmp.set(0,1,0).set(1,1,2).set(2,1,4);
209         matrix sol(A.solve(X, B));
210         for (unsigned ro=0; ro<3; ++ro)
211                 for (unsigned co=0; co<2; ++co)
212                         if (cmp(ro,co) != sol(ro,co))
213                                 result = 1;
214         if (result) {
215                 clog << "Solving " << A << " * " << X << " == " << B << endl
216                      << "erroneously returned " << sol << endl;
217         }
218         
219         return result;
220 }
221
222 static unsigned matrix_evalm()
223 {
224         unsigned result = 0;
225
226         matrix S(2, 2, lst(
227                 1, 2,
228                 3, 4
229         )), T(2, 2, lst(
230                 1, 1,
231                 2, -1
232         )), R(2, 2, lst(
233                 27, 14,
234                 36, 26
235         ));
236
237         ex e = ((S + T) * (S + 2*T));
238         ex f = e.evalm();
239         if (!f.is_equal(R)) {
240                 clog << "Evaluating " << e << " erroneously returned " << f << " instead of " << R << endl;
241                 result++;
242         }
243
244         return result;
245 }
246
247 static unsigned matrix_rank()
248 {
249         unsigned result = 0;
250         symbol x("x"), y("y");
251         matrix m(3,3);
252
253         // the zero matrix always has rank 0
254         if (m.rank() != 0) {
255                 clog << "The rank of " << m << " was not computed correctly." << endl;
256                 ++result;
257         }
258
259         // a trivial rank one example
260         m = 1, 0, 0,
261             2, 0, 0,
262             3, 0, 0;
263         if (m.rank() != 1) {
264                 clog << "The rank of " << m << " was not computed correctly." << endl;
265                 ++result;
266         }
267
268         // an example from Maple's help with rank two
269         m = x,  1,  0,
270             0,  0,  1,
271            x*y, y,  1;
272         if (m.rank() != 2) {
273                 clog << "The rank of " << m << " was not computed correctly." << endl;
274                 ++result;
275         }
276
277         // the 3x3 unit matrix has rank 3
278         m = ex_to<matrix>(unit_matrix(3,3));
279         if (m.rank() != 3) {
280                 clog << "The rank of " << m << " was not computed correctly." << endl;
281                 ++result;
282         }
283
284         return result;  
285 }
286
287 static unsigned matrix_misc()
288 {
289         unsigned result = 0;
290         matrix m1(2,2);
291         symbol a("a"), b("b"), c("c"), d("d"), e("e"), f("f");
292         m1.set(0,0,a).set(0,1,b);
293         m1.set(1,0,c).set(1,1,d);
294         ex tr = trace(m1);
295         
296         // check a simple trace
297         if (tr.compare(a+d)) {
298                 clog << "trace of 2x2 matrix " << m1
299                      << " erroneously returned " << tr << endl;
300                 ++result;
301         }
302         
303         // and two simple transpositions
304         matrix m2 = transpose(m1);
305         if (m2(0,0) != a || m2(0,1) != c || m2(1,0) != b || m2(1,1) != d) {
306                 clog << "transpose of 2x2 matrix " << m1
307                          << " erroneously returned " << m2 << endl;
308                 ++result;
309         }
310         matrix m3(3,2);
311         m3.set(0,0,a).set(0,1,b);
312         m3.set(1,0,c).set(1,1,d);
313         m3.set(2,0,e).set(2,1,f);
314         if (transpose(transpose(m3)) != m3) {
315                 clog << "transposing 3x2 matrix " << m3 << " twice"
316                      << " erroneously returned " << transpose(transpose(m3)) << endl;
317                 ++result;
318         }
319         
320         // produce a runtime-error by inverting a singular matrix and catch it
321         matrix m4(2,2);
322         matrix m5;
323         bool caught = false;
324         try {
325                 m5 = inverse(m4);
326         } catch (std::runtime_error err) {
327                 caught = true;
328         }
329         if (!caught) {
330                 cerr << "singular 2x2 matrix " << m4
331                      << " erroneously inverted to " << m5 << endl;
332                 ++result;
333         }
334         
335         return result;
336 }
337
338 unsigned exam_matrices()
339 {
340         unsigned result = 0;
341         
342         cout << "examining symbolic matrix manipulations" << flush;
343         
344         result += matrix_determinants();  cout << '.' << flush;
345         result += matrix_invert1();  cout << '.' << flush;
346         result += matrix_invert2();  cout << '.' << flush;
347         result += matrix_invert3();  cout << '.' << flush;
348         result += matrix_solve2();  cout << '.' << flush;
349         result += matrix_evalm();  cout << "." << flush;
350         result += matrix_rank();  cout << "." << flush;
351         result += matrix_misc();  cout << '.' << flush;
352         
353         return result;
354 }
355
356 int main(int argc, char** argv)
357 {
358         return exam_matrices();
359 }