Any complex number can be (un)archived properly.
[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-2008 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 <iostream>
25 #include "ginac.h"
26 using namespace std;
27 using namespace GiNaC;
28
29 static unsigned exam_powerlaws1()
30 {
31         // (x^a)^b = x^(a*b)
32         
33         symbol x("x");
34         symbol a("a");
35         symbol b("b");
36         
37         ex e1 = power(power(x,a), b);
38         if (!(is_exactly_a<power>(e1) &&
39               is_exactly_a<power>(e1.op(0)) &&
40               is_exactly_a<symbol>(e1.op(0).op(0)) &&
41               is_exactly_a<symbol>(e1.op(0).op(1)) &&
42               is_exactly_a<symbol>(e1.op(1)) &&
43               e1.is_equal(power(power(x,a),b)) )) {
44                 clog << "(x^a)^b, x,a,b symbolic wrong" << endl;
45                 clog << "returned: " << e1 << endl;
46                 return 1;
47         }
48         
49         ex e2 = e1.subs(a==1);
50         if (!(is_exactly_a<power>(e2) &&
51               is_exactly_a<symbol>(e2.op(0)) &&
52               is_exactly_a<symbol>(e2.op(1)) &&
53               e2.is_equal(power(x,b)) )) {
54                 clog << "(x^a)^b, x,b symbolic, a==1 wrong" << endl;
55                 clog << "returned: " << e2 << endl;
56                 return 1;
57         }
58         
59         ex e3 = e1.subs(a==-1);
60         if (!(is_exactly_a<power>(e3) &&
61               is_exactly_a<power>(e3.op(0)) &&
62               is_exactly_a<symbol>(e3.op(0).op(0)) &&
63               is_exactly_a<numeric>(e3.op(0).op(1)) &&
64               is_exactly_a<symbol>(e3.op(1)) &&
65               e3.is_equal(power(power(x,-1),b)) )) {
66                 clog << "(x^a)^b, x,b symbolic, a==-1 wrong" << endl;
67                 clog << "returned: " << e3 << endl;
68                 return 1;
69         }
70         
71         ex e4 = e1.subs(lst(a==-1, b==2.5));
72         if (!(is_exactly_a<power>(e4) &&
73               is_exactly_a<power>(e4.op(0)) &&
74               is_exactly_a<symbol>(e4.op(0).op(0)) &&
75               is_exactly_a<numeric>(e4.op(0).op(1)) &&
76               is_exactly_a<numeric>(e4.op(1)) &&
77               e4.is_equal(power(power(x,-1),2.5)) )) {
78                 clog << "(x^a)^b, x symbolic, a==-1, b==2.5 wrong" << endl;
79                 clog << "returned: " << e4 << endl;
80                 return 1;
81         }
82         
83         ex e5 = e1.subs(lst(a==-0.9, b==2.5));
84         if (!(is_exactly_a<power>(e5) &&
85               is_exactly_a<symbol>(e5.op(0)) &&
86               is_exactly_a<numeric>(e5.op(1)) &&
87               e5.is_equal(power(x,numeric(-0.9)*numeric(2.5))) )) {
88                 clog << "(x^a)^b, x symbolic, a==-0.9, b==2.5 wrong" << endl;
89                 clog << "returned: " << e5 << endl;
90                 return 1;
91         }
92         
93         ex e6 = e1.subs(lst(a==numeric(3)+numeric(5.3)*I, b==-5));
94         if (!(is_exactly_a<power>(e6) &&
95               is_exactly_a<symbol>(e6.op(0)) &&
96               is_exactly_a<numeric>(e6.op(1)) &&
97               e6.is_equal(power(x,numeric(-15)+numeric(5.3)*numeric(-5)*I)) )) {
98                 clog << "(x^a)^b, x symbolic, a==3+5.3*I, b==-5 wrong" << endl;
99                 clog << "returned: " << e6 << endl;
100                 return 1;
101         }
102         
103         return 0;
104 }
105
106 static unsigned exam_powerlaws2()
107 {
108         // (a*x)^b = a^b * x^b
109         
110         symbol x("x");
111         symbol a("a");
112         symbol b("b");
113         
114         ex e1 = power(a*x,b);
115         if (!(is_exactly_a<power>(e1) &&
116               is_exactly_a<mul>(e1.op(0)) &&
117               (e1.op(0).nops()==2) &&
118               is_exactly_a<symbol>(e1.op(0).op(0)) &&
119               is_exactly_a<symbol>(e1.op(0).op(1)) &&
120               is_exactly_a<symbol>(e1.op(1)) &&
121               e1.is_equal(power(a*x,b)) )) {
122                 clog << "(a*x)^b, x,a,b symbolic wrong" << endl;
123                 clog << "returned: " << e1 << endl;
124                 return 1;
125         }
126         
127         ex e2 = e1.subs(a==3);
128         if (!(is_exactly_a<power>(e2) &&
129               is_exactly_a<mul>(e2.op(0)) &&
130               (e2.op(0).nops()==2) &&
131               is_exactly_a<symbol>(e2.op(0).op(0)) &&
132               is_exactly_a<numeric>(e2.op(0).op(1)) &&
133               is_exactly_a<symbol>(e2.op(1)) &&
134               e2.is_equal(power(3*x,b)) )) {
135                 clog << "(a*x)^b, x,b symbolic, a==3 wrong" << endl;
136                 clog << "returned: " << e2 << endl;
137                 return 1;
138         }
139         
140         ex e3 = e1.subs(b==-3);
141         if (!(is_exactly_a<mul>(e3) &&
142               (e3.nops()==2) &&
143               is_exactly_a<power>(e3.op(0)) &&
144               is_exactly_a<power>(e3.op(1)) &&
145               e3.is_equal(power(a,-3)*power(x,-3)) )) {
146                 clog << "(a*x)^b, x,a symbolic, b==-3 wrong" << endl;
147                 clog << "returned: " << e3 << endl;
148                 return 1;
149         }
150         
151         ex e4 = e1.subs(b==4.5);
152         if (!(is_exactly_a<power>(e4) &&
153               is_exactly_a<mul>(e4.op(0)) &&
154               (e4.op(0).nops()==2) &&
155               is_exactly_a<symbol>(e4.op(0).op(0)) &&
156               is_exactly_a<symbol>(e4.op(0).op(1)) &&
157               is_exactly_a<numeric>(e4.op(1)) &&
158               e4.is_equal(power(a*x,4.5)) )) {
159                 clog << "(a*x)^b, x,a symbolic, b==4.5 wrong" << endl;
160                 clog << "returned: " << e4 << endl;
161                 return 1;
162         }
163         
164         ex e5 = e1.subs(lst(a==3.2, b==3+numeric(5)*I));
165         if (!(is_exactly_a<mul>(e5) &&
166               (e5.nops()==2) &&
167               is_exactly_a<power>(e5.op(0)) &&
168               is_exactly_a<numeric>(e5.op(1)) &&
169               e5.is_equal(power(x,3+numeric(5)*I)*
170                                           power(numeric(3.2),3+numeric(5)*I)) )) {
171                 clog << "(a*x)^b, x symbolic, a==3.2, b==3+5*I wrong" << endl;
172                 clog << "returned: " << e5 << endl;
173                 return 1;
174         }
175         
176         ex e6 = e1.subs(lst(a==-3.2, b==3+numeric(5)*I));
177         if (!(is_exactly_a<mul>(e6) &&
178               (e6.nops()==2) &&
179               is_exactly_a<power>(e6.op(0)) &&
180               is_exactly_a<numeric>(e6.op(1)) &&
181               e6.is_equal(power(-x,3+numeric(5)*I)*
182                                           power(numeric(3.2),3+numeric(5)*I)) )) {
183                 clog << "(a*x)^b, x symbolic, a==-3.2, b==3+5*I wrong" << endl;
184                 clog << "returned: " << e6 << endl;
185                 return 1;
186         }
187         
188         ex e7 = e1.subs(lst(a==3+numeric(5)*I, b==3.2));
189         if (!(is_exactly_a<power>(e7) &&
190               is_exactly_a<mul>(e7.op(0)) &&
191               (e7.op(0).nops()==2) &&
192               is_exactly_a<symbol>(e7.op(0).op(0)) &&
193               is_exactly_a<numeric>(e7.op(0).op(1)) &&
194               is_exactly_a<numeric>(e7.op(1)) &&
195               e7.is_equal(power((3+numeric(5)*I)*x,3.2)) )) {
196                 clog << "(a*x)^b, x symbolic, a==3+5*I, b==3.2 wrong" << endl;
197                 clog << "returned: " << e7 << endl;
198                 return 1;
199         }
200         
201         return 0;
202 }
203
204 static unsigned exam_powerlaws3()
205 {
206         // numeric evaluation
207
208         ex e1 = power(numeric(4),numeric(1,2));
209         if (e1 != 2) {
210                 clog << "4^(1/2) wrongly returned " << e1 << endl;
211                 return 1;
212         }
213         
214         ex e2 = power(numeric(27),numeric(2,3));
215         if (e2 != 9) {
216                 clog << "27^(2/3) wrongly returned " << e2 << endl;
217                 return 1;
218         }
219         
220         ex e3 = power(numeric(5),numeric(1,2));
221         if (!(is_exactly_a<power>(e3) &&
222               e3.op(0).is_equal(numeric(5)) &&
223               e3.op(1).is_equal(numeric(1,2)))) {
224                 clog << "5^(1/2) wrongly returned " << e3 << endl;
225                 return 1;
226         }
227         
228         ex e4 = power(numeric(5),evalf(numeric(1,2)));
229         if (!(is_exactly_a<numeric>(e4))) {
230                 clog << "5^(0.5) wrongly returned " << e4 << endl;
231                 return 1;
232         }
233         
234         ex e5 = power(evalf(numeric(5)),numeric(1,2));
235         if (!(is_exactly_a<numeric>(e5))) {
236                 clog << "5.0^(1/2) wrongly returned " << e5 << endl;
237                 return 1;
238         }
239         
240         return 0;
241 }
242
243 static unsigned exam_powerlaws4()
244 {
245         // test for mul::eval()
246         
247         symbol a("a");
248         symbol b("b");
249         symbol c("c");
250         
251         ex f1 = power(a*b,ex(1)/ex(2));
252         ex f2 = power(a*b,ex(3)/ex(2));
253         ex f3 = c;
254         
255         exvector v;
256         v.push_back(f1);
257         v.push_back(f2);
258         v.push_back(f3);
259         ex e1 = mul(v);
260         if (e1!=a*a*b*b*c) {
261                 clog << "(a*b)^(1/2)*(a*b)^(3/2)*c wrongly returned " << e1 << endl;
262                 return 1;
263         }
264         
265         return 0;
266 }
267
268 static unsigned exam_powerlaws5()
269 {
270         // cabinet of slightly pathological cases
271         
272         symbol a("a");
273         
274         ex e1 = pow(1,a);
275         if (e1 != 1) {
276                 clog << "1^a wrongly returned " << e1 << endl;
277                 return 1;
278         }
279         
280         ex e2 = pow(0,a);
281         if (!(is_exactly_a<power>(e2))) {
282                 clog << "0^a was evaluated to " << e2
283                      << " though nothing is known about a." << endl;
284                 return 1;
285         }
286         
287         return 0;
288 }
289
290 unsigned exam_powerlaws()
291 {
292         unsigned result = 0;
293         
294         cout << "examining power laws" << flush;
295         
296         result += exam_powerlaws1();  cout << '.' << flush;
297         result += exam_powerlaws2();  cout << '.' << flush;
298         result += exam_powerlaws3();  cout << '.' << flush;
299         result += exam_powerlaws4();  cout << '.' << flush;
300         result += exam_powerlaws5();  cout << '.' << flush;
301         
302         return result;
303 }
304
305 int main(int argc, char** argv)
306 {
307         return exam_powerlaws();
308 }