- enforced GiNaC coding standards :-)
[ginac.git] / check / powerlaws.cpp
1 /** @file 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  *  GiNaC Copyright (C) 1999 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
21  */
22
23 #include <ginac/ginac.h>
24
25 static unsigned powerlaws1(void)
26 {
27     // (x^a)^b = x^(a*b)
28     
29     symbol x("x");
30     symbol a("a");
31     symbol b("b");
32     
33     ex e1=power(power(x,a),b);
34     if (!(is_ex_exactly_of_type(e1,power) &&
35           is_ex_exactly_of_type(e1.op(0),power) &&
36           is_ex_exactly_of_type(e1.op(0).op(0),symbol) &&
37           is_ex_exactly_of_type(e1.op(0).op(1),symbol) &&
38           is_ex_exactly_of_type(e1.op(1),symbol) &&
39           e1.is_equal(power(power(x,a),b)) )) {
40         clog << "(x^a)^b, x,a,b symbolic wrong" << endl;
41         clog << "returned: " << e1 << endl;
42         return 1;
43     }
44     
45     ex e2=e1.subs(a==1);
46     if (!(is_ex_exactly_of_type(e2,power) &&
47           is_ex_exactly_of_type(e2.op(0),symbol) &&
48           is_ex_exactly_of_type(e2.op(1),symbol) &&
49           e2.is_equal(power(x,b)) )) {
50         clog << "(x^a)^b, x,b symbolic, a==1 wrong" << endl;
51         clog << "returned: " << e2 << endl;
52         return 1;
53     }
54     
55     ex e3=e1.subs(a==-1);
56     if (!(is_ex_exactly_of_type(e3,power) &&
57           is_ex_exactly_of_type(e3.op(0),power) &&
58           is_ex_exactly_of_type(e3.op(0).op(0),symbol) &&
59           is_ex_exactly_of_type(e3.op(0).op(1),numeric) &&
60           is_ex_exactly_of_type(e3.op(1),symbol) &&
61           e3.is_equal(power(power(x,-1),b)) )) {
62         clog << "(x^a)^b, x,b symbolic, a==-1 wrong" << endl;
63         clog << "returned: " << e3 << endl;
64         return 1;
65     }
66     
67     ex e4=e1.subs(lst(a==-1,b==2.5));
68     if (!(is_ex_exactly_of_type(e4,power) &&
69           is_ex_exactly_of_type(e4.op(0),power) &&
70           is_ex_exactly_of_type(e4.op(0).op(0),symbol) &&
71           is_ex_exactly_of_type(e4.op(0).op(1),numeric) &&
72           is_ex_exactly_of_type(e4.op(1),numeric) &&
73           e4.is_equal(power(power(x,-1),2.5)) )) {
74         clog << "(x^a)^b, x symbolic, a==-1, b==2.5 wrong" << endl;
75         clog << "returned: " << e4 << endl;
76         return 1;
77     }
78     
79     ex e5=e1.subs(lst(a==-0.9,b==2.5));
80     if (!(is_ex_exactly_of_type(e5,power) &&
81           is_ex_exactly_of_type(e5.op(0),symbol) &&
82           is_ex_exactly_of_type(e5.op(1),numeric) &&
83           e5.is_equal(power(x,numeric(-0.9)*numeric(2.5))) )) {
84         clog << "(x^a)^b, x symbolic, a==-0.9, b==2.5 wrong" << endl;
85         clog << "returned: " << e5 << endl;
86         return 1;
87     }
88     
89     ex e6=e1.subs(lst(a==numeric(3)+numeric(5.3)*I,b==-5));
90     if (!(is_ex_exactly_of_type(e6,power) &&
91           is_ex_exactly_of_type(e6.op(0),symbol) &&
92           is_ex_exactly_of_type(e6.op(1),numeric) &&
93           e6.is_equal(power(x,numeric(-15)+numeric(5.3)*numeric(-5)*I)) )) {
94         clog << "(x^a)^b, x symbolic, a==3+5.3*I, b==-5 wrong" << endl;
95         clog << "returned: " << e6 << endl;
96         return 1;
97     }
98     return 0;
99 }
100
101 static unsigned powerlaws2(void)
102 {
103     // (a*x)^b = a^b * x^b
104     
105     symbol x("x");
106     symbol a("a");
107     symbol b("b");
108     
109     ex e1=power(a*x,b);
110     if (!(is_ex_exactly_of_type(e1,power) &&
111           is_ex_exactly_of_type(e1.op(0),mul) &&
112           (e1.op(0).nops()==2) &&
113           is_ex_exactly_of_type(e1.op(0).op(0),symbol) &&
114           is_ex_exactly_of_type(e1.op(0).op(1),symbol) &&
115           is_ex_exactly_of_type(e1.op(1),symbol) &&
116           e1.is_equal(power(a*x,b)) )) {
117         clog << "(a*x)^b, x,a,b symbolic wrong" << endl;
118         clog << "returned: " << e1 << endl;
119         return 1;
120     }
121     
122     ex e2=e1.subs(a==3);
123     if (!(is_ex_exactly_of_type(e2,power) &&
124           is_ex_exactly_of_type(e2.op(0),mul) &&
125           (e2.op(0).nops()==2) &&
126           is_ex_exactly_of_type(e2.op(0).op(0),symbol) &&
127           is_ex_exactly_of_type(e2.op(0).op(1),numeric) &&
128           is_ex_exactly_of_type(e2.op(1),symbol) &&
129           e2.is_equal(power(3*x,b)) )) {
130         clog << "(a*x)^b, x,b symbolic, a==3 wrong" << endl;
131         clog << "returned: " << e2 << endl;
132         return 1;
133     }
134     
135     ex e3=e1.subs(b==-3);
136     if (!(is_ex_exactly_of_type(e3,mul) &&
137           (e3.nops()==2) &&
138           is_ex_exactly_of_type(e3.op(0),power) &&
139           is_ex_exactly_of_type(e3.op(1),power) &&
140           e3.is_equal(power(a,-3)*power(x,-3)) )) {
141         clog << "(a*x)^b, x,a symbolic, b==-3 wrong" << endl;
142         clog << "returned: " << e3 << endl;
143         return 1;
144     }
145     
146     ex e4=e1.subs(b==4.5);
147     if (!(is_ex_exactly_of_type(e4,power) &&
148           is_ex_exactly_of_type(e4.op(0),mul) &&
149           (e4.op(0).nops()==2) &&
150           is_ex_exactly_of_type(e4.op(0).op(0),symbol) &&
151           is_ex_exactly_of_type(e4.op(0).op(1),symbol) &&
152           is_ex_exactly_of_type(e4.op(1),numeric) &&
153           e4.is_equal(power(a*x,4.5)) )) {
154         clog << "(a*x)^b, x,a symbolic, b==4.5 wrong" << endl;
155         clog << "returned: " << e4 << endl;
156         return 1;
157     }
158     
159     ex e5=e1.subs(lst(a==3.2,b==3+numeric(5)*I));
160     if (!(is_ex_exactly_of_type(e5,mul) &&
161           (e5.nops()==2) &&
162           is_ex_exactly_of_type(e5.op(0),power) &&
163           is_ex_exactly_of_type(e5.op(1),numeric) &&
164           e5.is_equal(power(x,3+numeric(5)*I)*
165                       power(numeric(3.2),3+numeric(5)*I)) )) {
166         clog << "(a*x)^b, x symbolic, a==3.2, b==3+5*I wrong" << endl;
167         clog << "returned: " << e5 << endl;
168         return 1;
169     }
170     
171     ex e6=e1.subs(lst(a==-3.2,b==3+numeric(5)*I));
172     if (!(is_ex_exactly_of_type(e6,mul) &&
173           (e6.nops()==2) &&
174           is_ex_exactly_of_type(e6.op(0),power) &&
175           is_ex_exactly_of_type(e6.op(1),numeric) &&
176           e6.is_equal(power(-x,3+numeric(5)*I)*
177                       power(numeric(3.2),3+numeric(5)*I)) )) {
178         clog << "(a*x)^b, x symbolic, a==-3.2, b==3+5*I wrong" << endl;
179         clog << "returned: " << e6 << endl;
180         return 1;
181     }
182     
183     ex e7=e1.subs(lst(a==3+numeric(5)*I,b==3.2));
184     if (!(is_ex_exactly_of_type(e7,power) &&
185           is_ex_exactly_of_type(e7.op(0),mul) &&
186           (e7.op(0).nops()==2) &&
187           is_ex_exactly_of_type(e7.op(0).op(0),symbol) &&
188           is_ex_exactly_of_type(e7.op(0).op(1),numeric) &&
189           is_ex_exactly_of_type(e7.op(1),numeric) &&
190           e7.is_equal(power((3+numeric(5)*I)*x,3.2)) )) {
191         clog << "(a*x)^b, x symbolic, a==3+5*I, b==3.2 wrong" << endl;
192         clog << "returned: " << e7 << endl;
193         return 1;
194     }
195     
196     return 0;
197 }
198
199 static unsigned powerlaws3(void)
200 {
201     // numeric evaluation
202
203     ex e1=power(numeric(4),numeric(1)/numeric(2));
204     if (e1 != 2) {
205         clog << "4^(1/2) wrongly returned " << e1 << endl;
206         return 1;
207     }
208     
209     ex e2=power(numeric(27),numeric(2)/numeric(3));
210     if (e2 != 9) {
211         clog << "27^(2/3) wrongly returned " << e2 << endl;
212         return 1;
213     }
214
215     ex e3=power(numeric(5),numeric(1)/numeric(2));
216     if (!(is_ex_exactly_of_type(e3,power) &&
217           e3.op(0).is_equal(numeric(5)) &&
218           e3.op(1).is_equal(numeric(1)/numeric(2)) )) {
219         clog << "5^(1/2) wrongly returned " << e3 << endl;
220         return 1;
221     }
222     
223     ex e4=power(numeric(5),evalf(numeric(1)/numeric(2)));
224     if (!(is_ex_exactly_of_type(e4,numeric))) {
225         clog << "5^(0.5) wrongly returned " << e4 << endl;
226         return 1;
227     }
228     
229     ex e5=power(evalf(numeric(5)),numeric(1)/numeric(2));
230     if (!(is_ex_exactly_of_type(e5,numeric))) {
231         clog << "5.0^(1/2) wrongly returned " << e5 << endl;
232         return 1;
233     }
234     
235     return 0;
236 }
237
238 static unsigned powerlaws4(void)
239 {
240     // test for mul::eval()
241
242     symbol a("a");
243     symbol b("b");
244     symbol c("c");
245     
246     ex f1=power(a*b,ex(1)/ex(2));
247     ex f2=power(a*b,ex(3)/ex(2));
248     ex f3=c;
249
250     exvector v;
251     v.push_back(f1);
252     v.push_back(f2);
253     v.push_back(f3);
254     ex e1=mul(v);
255     if (e1!=a*a*b*b*c) {
256         clog << "(a*b)^(1/2)*(a*b)^(3/2)*c wrongly returned " << e1 << endl;
257         return 1;
258     }
259     return 0;
260 }
261
262 unsigned powerlaws(void)
263 {
264     unsigned result = 0;
265     
266     cout << "checking power laws..." << flush;
267     clog << "---------power laws:" << endl;
268     
269     result += powerlaws1();
270     result += powerlaws2();
271     result += powerlaws3();
272     result += powerlaws4();
273     
274     if (!result) {
275         cout << " passed ";
276         clog << "(no output)" << endl;
277     } else {
278         cout << " failed ";
279     }
280     return result;
281 }