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-2000 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 }