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 }