#ifndef around namespace GiNaC { }
[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 /*
7  *  GiNaC Copyright (C) 1999 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
22  */
23
24 #include <ginac/ginac.h>
25
26 #ifndef NO_GINAC_NAMESPACE
27 using namespace GiNaC;
28 #endif // ndef NO_GINAC_NAMESPACE
29
30 static unsigned powerlaws1(void)
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_ex_exactly_of_type(e1,power) &&
40           is_ex_exactly_of_type(e1.op(0),power) &&
41           is_ex_exactly_of_type(e1.op(0).op(0),symbol) &&
42           is_ex_exactly_of_type(e1.op(0).op(1),symbol) &&
43           is_ex_exactly_of_type(e1.op(1),symbol) &&
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_ex_exactly_of_type(e2,power) &&
52           is_ex_exactly_of_type(e2.op(0),symbol) &&
53           is_ex_exactly_of_type(e2.op(1),symbol) &&
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_ex_exactly_of_type(e3,power) &&
62           is_ex_exactly_of_type(e3.op(0),power) &&
63           is_ex_exactly_of_type(e3.op(0).op(0),symbol) &&
64           is_ex_exactly_of_type(e3.op(0).op(1),numeric) &&
65           is_ex_exactly_of_type(e3.op(1),symbol) &&
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_ex_exactly_of_type(e4,power) &&
74           is_ex_exactly_of_type(e4.op(0),power) &&
75           is_ex_exactly_of_type(e4.op(0).op(0),symbol) &&
76           is_ex_exactly_of_type(e4.op(0).op(1),numeric) &&
77           is_ex_exactly_of_type(e4.op(1),numeric) &&
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_ex_exactly_of_type(e5,power) &&
86           is_ex_exactly_of_type(e5.op(0),symbol) &&
87           is_ex_exactly_of_type(e5.op(1),numeric) &&
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_ex_exactly_of_type(e6,power) &&
96           is_ex_exactly_of_type(e6.op(0),symbol) &&
97           is_ex_exactly_of_type(e6.op(1),numeric) &&
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     return 0;
104 }
105
106 static unsigned powerlaws2(void)
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_ex_exactly_of_type(e1,power) &&
116           is_ex_exactly_of_type(e1.op(0),mul) &&
117           (e1.op(0).nops()==2) &&
118           is_ex_exactly_of_type(e1.op(0).op(0),symbol) &&
119           is_ex_exactly_of_type(e1.op(0).op(1),symbol) &&
120           is_ex_exactly_of_type(e1.op(1),symbol) &&
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_ex_exactly_of_type(e2,power) &&
129           is_ex_exactly_of_type(e2.op(0),mul) &&
130           (e2.op(0).nops()==2) &&
131           is_ex_exactly_of_type(e2.op(0).op(0),symbol) &&
132           is_ex_exactly_of_type(e2.op(0).op(1),numeric) &&
133           is_ex_exactly_of_type(e2.op(1),symbol) &&
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_ex_exactly_of_type(e3,mul) &&
142           (e3.nops()==2) &&
143           is_ex_exactly_of_type(e3.op(0),power) &&
144           is_ex_exactly_of_type(e3.op(1),power) &&
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_ex_exactly_of_type(e4,power) &&
153           is_ex_exactly_of_type(e4.op(0),mul) &&
154           (e4.op(0).nops()==2) &&
155           is_ex_exactly_of_type(e4.op(0).op(0),symbol) &&
156           is_ex_exactly_of_type(e4.op(0).op(1),symbol) &&
157           is_ex_exactly_of_type(e4.op(1),numeric) &&
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_ex_exactly_of_type(e5,mul) &&
166           (e5.nops()==2) &&
167           is_ex_exactly_of_type(e5.op(0),power) &&
168           is_ex_exactly_of_type(e5.op(1),numeric) &&
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_ex_exactly_of_type(e6,mul) &&
178           (e6.nops()==2) &&
179           is_ex_exactly_of_type(e6.op(0),power) &&
180           is_ex_exactly_of_type(e6.op(1),numeric) &&
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_ex_exactly_of_type(e7,power) &&
190           is_ex_exactly_of_type(e7.op(0),mul) &&
191           (e7.op(0).nops()==2) &&
192           is_ex_exactly_of_type(e7.op(0).op(0),symbol) &&
193           is_ex_exactly_of_type(e7.op(0).op(1),numeric) &&
194           is_ex_exactly_of_type(e7.op(1),numeric) &&
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 powerlaws3(void)
205 {
206     // numeric evaluation
207
208     ex e1=power(numeric(4),numeric(1)/numeric(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)/numeric(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)/numeric(2));
221     if (!(is_ex_exactly_of_type(e3,power) &&
222           e3.op(0).is_equal(numeric(5)) &&
223           e3.op(1).is_equal(numeric(1)/numeric(2)) )) {
224         clog << "5^(1/2) wrongly returned " << e3 << endl;
225         return 1;
226     }
227     
228     ex e4=power(numeric(5),evalf(numeric(1)/numeric(2)));
229     if (!(is_ex_exactly_of_type(e4,numeric))) {
230         clog << "5^(0.5) wrongly returned " << e4 << endl;
231         return 1;
232     }
233     
234     ex e5=power(evalf(numeric(5)),numeric(1)/numeric(2));
235     if (!(is_ex_exactly_of_type(e5,numeric))) {
236         clog << "5.0^(1/2) wrongly returned " << e5 << endl;
237         return 1;
238     }
239     
240     return 0;
241 }
242
243 static unsigned powerlaws4(void)
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     return 0;
265 }
266
267 unsigned powerlaws(void)
268 {
269     unsigned result = 0;
270     
271     cout << "checking power laws..." << flush;
272     clog << "---------power laws:" << endl;
273     
274     result += powerlaws1();
275     result += powerlaws2();
276     result += powerlaws3();
277     result += powerlaws4();
278     
279     if (!result) {
280         cout << " passed ";
281         clog << "(no output)" << endl;
282     } else {
283         cout << " failed ";
284     }
285     return result;
286 }