Happy New Year!
[ginac.git] / check / exam_powerlaws.cpp
1 /** @file exam_powerlaws.cpp
2  *
3  *  Tests for power laws.  You shouldn't try to draw much inspiration from
4  *  this code, it is a sanity check rather deeply rooted in GiNaC's classes. */
5
6 /*
7  *  GiNaC Copyright (C) 1999-2019 Johannes Gutenberg University Mainz, Germany
8  *
9  *  This program is free software; you can redistribute it and/or modify
10  *  it under the terms of the GNU General Public License as published by
11  *  the Free Software Foundation; either version 2 of the License, or
12  *  (at your option) any later version.
13  *
14  *  This program is distributed in the hope that it will be useful,
15  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17  *  GNU General Public License for more details.
18  *
19  *  You should have received a copy of the GNU General Public License
20  *  along with this program; if not, write to the Free Software
21  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
22  */
23
24 #include "ginac.h"
25 using namespace GiNaC;
26
27 #include <iostream>
28 using namespace std;
29
30 static unsigned exam_powerlaws1()
31 {
32         // (x^a)^b = x^(a*b)
33         
34         symbol x("x");
35         symbol a("a");
36         symbol b("b");
37         
38         ex e1 = power(power(x,a), b);
39         if (!(is_exactly_a<power>(e1) &&
40               is_exactly_a<power>(e1.op(0)) &&
41               is_exactly_a<symbol>(e1.op(0).op(0)) &&
42               is_exactly_a<symbol>(e1.op(0).op(1)) &&
43               is_exactly_a<symbol>(e1.op(1)) &&
44               e1.is_equal(power(power(x,a),b)) )) {
45                 clog << "(x^a)^b, x,a,b symbolic wrong" << endl;
46                 clog << "returned: " << e1 << endl;
47                 return 1;
48         }
49         
50         ex e2 = e1.subs(a==1);
51         if (!(is_exactly_a<power>(e2) &&
52               is_exactly_a<symbol>(e2.op(0)) &&
53               is_exactly_a<symbol>(e2.op(1)) &&
54               e2.is_equal(power(x,b)) )) {
55                 clog << "(x^a)^b, x,b symbolic, a==1 wrong" << endl;
56                 clog << "returned: " << e2 << endl;
57                 return 1;
58         }
59         
60         ex e3 = e1.subs(a==-1);
61         if (!(is_exactly_a<power>(e3) &&
62               is_exactly_a<power>(e3.op(0)) &&
63               is_exactly_a<symbol>(e3.op(0).op(0)) &&
64               is_exactly_a<numeric>(e3.op(0).op(1)) &&
65               is_exactly_a<symbol>(e3.op(1)) &&
66               e3.is_equal(power(power(x,-1),b)) )) {
67                 clog << "(x^a)^b, x,b symbolic, a==-1 wrong" << endl;
68                 clog << "returned: " << e3 << endl;
69                 return 1;
70         }
71         
72         ex e4 = e1.subs(lst{a==-1, b==-2.5});
73         if (!(is_exactly_a<power>(e4) &&
74               is_exactly_a<power>(e4.op(0)) &&
75               is_exactly_a<symbol>(e4.op(0).op(0)) &&
76               is_exactly_a<numeric>(e4.op(0).op(1)) &&
77               is_exactly_a<numeric>(e4.op(1)) &&
78               e4.is_equal(power(power(x,-1),-2.5)) )) {
79                 clog << "(x^a)^b, x symbolic, a==-1, b==-2.5 wrong" << endl;
80                 clog << "returned: " << e4 << endl;
81                 return 1;
82         }
83         
84         ex e5 = e1.subs(lst{a==-0.9, b==2.5});
85         if (!(is_exactly_a<power>(e5) &&
86               is_exactly_a<symbol>(e5.op(0)) &&
87               is_exactly_a<numeric>(e5.op(1)) &&
88               e5.is_equal(power(x,numeric(-0.9)*numeric(2.5))) )) {
89                 clog << "(x^a)^b, x symbolic, a==-0.9, b==2.5 wrong" << endl;
90                 clog << "returned: " << e5 << endl;
91                 return 1;
92         }
93         
94         ex e6 = e1.subs(lst{a==numeric(3)+numeric(5.3)*I, b==-5});
95         if (!(is_exactly_a<power>(e6) &&
96               is_exactly_a<symbol>(e6.op(0)) &&
97               is_exactly_a<numeric>(e6.op(1)) &&
98               e6.is_equal(power(x,numeric(-15)+numeric(5.3)*numeric(-5)*I)) )) {
99                 clog << "(x^a)^b, x symbolic, a==3+5.3*I, b==-5 wrong" << endl;
100                 clog << "returned: " << e6 << endl;
101                 return 1;
102         }
103         
104         return 0;
105 }
106
107 static unsigned exam_powerlaws2()
108 {
109         // (a*x)^b = a^b * x^b
110         
111         symbol x("x");
112         symbol a("a");
113         symbol b("b");
114         
115         ex e1 = power(a*x,b);
116         if (!(is_exactly_a<power>(e1) &&
117               is_exactly_a<mul>(e1.op(0)) &&
118               (e1.op(0).nops()==2) &&
119               is_exactly_a<symbol>(e1.op(0).op(0)) &&
120               is_exactly_a<symbol>(e1.op(0).op(1)) &&
121               is_exactly_a<symbol>(e1.op(1)) &&
122               e1.is_equal(power(a*x,b)) )) {
123                 clog << "(a*x)^b, x,a,b symbolic wrong" << endl;
124                 clog << "returned: " << e1 << endl;
125                 return 1;
126         }
127         
128         ex e2 = e1.subs(a==3);
129         if (!(is_exactly_a<power>(e2) &&
130               is_exactly_a<mul>(e2.op(0)) &&
131               (e2.op(0).nops()==2) &&
132               is_exactly_a<symbol>(e2.op(0).op(0)) &&
133               is_exactly_a<numeric>(e2.op(0).op(1)) &&
134               is_exactly_a<symbol>(e2.op(1)) &&
135               e2.is_equal(power(3*x,b)) )) {
136                 clog << "(a*x)^b, x,b symbolic, a==3 wrong" << endl;
137                 clog << "returned: " << e2 << endl;
138                 return 1;
139         }
140         
141         ex e3 = e1.subs(b==-3);
142         if (!(is_exactly_a<mul>(e3) &&
143               (e3.nops()==2) &&
144               is_exactly_a<power>(e3.op(0)) &&
145               is_exactly_a<power>(e3.op(1)) &&
146               e3.is_equal(power(a,-3)*power(x,-3)) )) {
147                 clog << "(a*x)^b, x,a symbolic, b==-3 wrong" << endl;
148                 clog << "returned: " << e3 << endl;
149                 return 1;
150         }
151         
152         ex e4 = e1.subs(b==4.5);
153         if (!(is_exactly_a<power>(e4) &&
154               is_exactly_a<mul>(e4.op(0)) &&
155               (e4.op(0).nops()==2) &&
156               is_exactly_a<symbol>(e4.op(0).op(0)) &&
157               is_exactly_a<symbol>(e4.op(0).op(1)) &&
158               is_exactly_a<numeric>(e4.op(1)) &&
159               e4.is_equal(power(a*x,4.5)) )) {
160                 clog << "(a*x)^b, x,a symbolic, b==4.5 wrong" << endl;
161                 clog << "returned: " << e4 << endl;
162                 return 1;
163         }
164         
165         ex e5 = e1.subs(lst{a==3.2, b==3+numeric(5)*I});
166         if (!(is_exactly_a<mul>(e5) &&
167               (e5.nops()==2) &&
168               is_exactly_a<power>(e5.op(0)) &&
169               is_exactly_a<numeric>(e5.op(1)) &&
170               e5.is_equal(power(x,3+numeric(5)*I)*
171                                           power(numeric(3.2),3+numeric(5)*I)) )) {
172                 clog << "(a*x)^b, x symbolic, a==3.2, b==3+5*I wrong" << endl;
173                 clog << "returned: " << e5 << endl;
174                 return 1;
175         }
176         
177         ex e6 = e1.subs(lst{a==-3.2, b==3+numeric(5)*I});
178         if (!(is_exactly_a<mul>(e6) &&
179               (e6.nops()==2) &&
180               is_exactly_a<power>(e6.op(0)) &&
181               is_exactly_a<numeric>(e6.op(1)) &&
182               e6.is_equal(power(-x,3+numeric(5)*I)*
183                                           power(numeric(3.2),3+numeric(5)*I)) )) {
184                 clog << "(a*x)^b, x symbolic, a==-3.2, b==3+5*I wrong" << endl;
185                 clog << "returned: " << e6 << endl;
186                 return 1;
187         }
188         
189         ex e7 = e1.subs(lst{a==3+numeric(5)*I, b==3.2});
190         if (!(is_exactly_a<power>(e7) &&
191               is_exactly_a<mul>(e7.op(0)) &&
192               (e7.op(0).nops()==2) &&
193               is_exactly_a<symbol>(e7.op(0).op(0)) &&
194               is_exactly_a<numeric>(e7.op(0).op(1)) &&
195               is_exactly_a<numeric>(e7.op(1)) &&
196               e7.is_equal(power((3+numeric(5)*I)*x,3.2)) )) {
197                 clog << "(a*x)^b, x symbolic, a==3+5*I, b==3.2 wrong" << endl;
198                 clog << "returned: " << e7 << endl;
199                 return 1;
200         }
201         
202         return 0;
203 }
204
205 static unsigned exam_powerlaws3()
206 {
207         // numeric evaluation
208
209         ex e1 = power(numeric(4),numeric(1,2));
210         if (e1 != 2) {
211                 clog << "4^(1/2) wrongly returned " << e1 << endl;
212                 return 1;
213         }
214         
215         ex e2 = power(numeric(27),numeric(2,3));
216         if (e2 != 9) {
217                 clog << "27^(2/3) wrongly returned " << e2 << endl;
218                 return 1;
219         }
220         
221         ex e3 = power(numeric(5),numeric(1,2));
222         if (!(is_exactly_a<power>(e3) &&
223               e3.op(0).is_equal(numeric(5)) &&
224               e3.op(1).is_equal(numeric(1,2)))) {
225                 clog << "5^(1/2) wrongly returned " << e3 << endl;
226                 return 1;
227         }
228         
229         ex e4 = power(numeric(5),evalf(numeric(1,2)));
230         if (!(is_exactly_a<numeric>(e4))) {
231                 clog << "5^(0.5) wrongly returned " << e4 << endl;
232                 return 1;
233         }
234         
235         ex e5 = power(evalf(numeric(5)),numeric(1,2));
236         if (!(is_exactly_a<numeric>(e5))) {
237                 clog << "5.0^(1/2) wrongly returned " << e5 << endl;
238                 return 1;
239         }
240         
241         return 0;
242 }
243
244 static unsigned exam_powerlaws4()
245 {
246         // test for mul::eval()
247         
248         symbol a("a");
249         symbol b("b");
250         symbol c("c");
251         
252         ex f1 = power(a*b,ex(1)/ex(2));
253         ex f2 = power(a*b,ex(3)/ex(2));
254         ex f3 = c;
255         
256         exvector v;
257         v.push_back(f1);
258         v.push_back(f2);
259         v.push_back(f3);
260         ex e1 = mul(v);
261         if (e1!=a*a*b*b*c) {
262                 clog << "(a*b)^(1/2)*(a*b)^(3/2)*c wrongly returned " << e1 << endl;
263                 return 1;
264         }
265         
266         return 0;
267 }
268
269 static unsigned exam_powerlaws5()
270 {
271         // cabinet of slightly pathological cases
272         
273         symbol a("a");
274         
275         ex e1 = pow(1,a);
276         if (e1 != 1) {
277                 clog << "1^a wrongly returned " << e1 << endl;
278                 return 1;
279         }
280         
281         ex e2 = pow(0,a);
282         if (!(is_exactly_a<power>(e2))) {
283                 clog << "0^a was evaluated to " << e2
284                      << " though nothing is known about a." << endl;
285                 return 1;
286         }
287         
288         return 0;
289 }
290
291 static unsigned exam_powerlaws6()
292 {
293         // check expansion rules for positive symbols
294
295         symbol a("a");
296         symbol b("b");
297         symbol c("c");
298         realsymbol x("x");
299         realsymbol y("y");
300         possymbol p("p");
301         possymbol q("q");
302         numeric half=numeric(1,2);
303
304         ex e1 = pow(5*pow(3*a*b*x*y*p*q,2),7*half*c).expand();
305         ex e2 = pow(p,7*c)*pow(q,7*c)*pow(pow(a*b*x*y,2),numeric(7,2)*c)*pow(45,numeric(7,2)*c);
306         if (!e1.is_equal(e2)) {
307                 clog << "Could not expand exponents with positive bases in " << e1 << endl;
308                 return 1;
309         }
310
311         ex e3 = pow(-pow(-a*x*p,3)*pow(b*y*p,3),half*c).expand().normal();
312         ex e4 = pow(p,3*c)*pow(pow(a*b*x*y,3),half*c);
313
314         if (!e3.is_equal(e4)) {
315                 clog << "Could not expand exponents with positive bases in " << e3 << endl;
316                 return 1;
317         }
318
319         return 0;
320 }
321
322 unsigned exam_powerlaws()
323 {
324         unsigned result = 0;
325         
326         cout << "examining power laws" << flush;
327         
328         result += exam_powerlaws1();  cout << '.' << flush;
329         result += exam_powerlaws2();  cout << '.' << flush;
330         result += exam_powerlaws3();  cout << '.' << flush;
331         result += exam_powerlaws4();  cout << '.' << flush;
332         result += exam_powerlaws5();  cout << '.' << flush;
333         result += exam_powerlaws6();  cout << '.' << flush;
334         
335         return result;
336 }
337
338 int main(int argc, char** argv)
339 {
340         return exam_powerlaws();
341 }