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 }